NAVAL  POSTGRADUATE  SCHOOL 
MONTEREY,  CALIFORNIA 


THESIS 


DETECTION  AND  IDENTIFICATION 

OF 
CYCLOSTATIONARY  SIGNALS 


by 


Evandro  Luiz  da  Costa 


March,  1996 


Thesis  Co-Advisors: 


Ralph  Hippenstiel 
Roberto  Cristi 


Approved  for  public  release;  distribution  is  unlimited. 


DUDLEY  KNOX  LIBRARY 

-GRADUATE  SCHOOL 


REPORT  DOCUMENTATION  PAGE 


Form  Approved  OMB  No.  0704-0188 


Public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instruction,  searching 
existing  data  sources,  gathering  and  maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this 
burden  estimate  or  any  other  aspect  of  this  collection  of  information,  including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services, 
Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington,  VA  22202-4302,  and  to  the  Office  of  Management  and 
Budget,  Paperwork  Reduction  Project  (0704-0188)  Washington  DC  20503. 


1 .      AGENCY  USE  ONLY  (Leave  blank) 


2.      REPORT  DATE 

March  1996 


REPORT  TYPE  AND  DATES  COVERED 

Master's  Thesis 


TITLE  AND  SUBTITLE 

DETECTION  AND  IDENTIFICATION  OF  CYCLOSTATIONARY 
SIGNALS 


AUTHOR(S) 

Da  Costa,  Evandro  L. 


FUNDING  NUMBERS 


7.      PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Naval  Postgraduate  School 
Monterey  CA  93943-5000 


PERFORMING 
ORGANIZATION 
REPORT  NUMBER 


9.      SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 


10.    SPONSORING/MONITORING 
AGENCY  REPORT  NUMBER 


11.  SUPPLEMENTARY  NOTES 

The  views  expressed  in  this  thesis  are  those  of  the  author  and  do  not  reflect  the  official  policy  or 
position  of  the  Department  of  Defense  or  the  U.S.  Government. 


12a.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  is  unlimited. 


12b.  DISTRIBUTION  CODE 


13.    ABSTRACT  (maximum  200  words) 

Propeller  noise  can  be  modeled  as  an  amplitude  modulated  (AM)  signal.  Cyclic  Spectral 
Analysis  has  been  used  successfully  to  detect  the  presence  of  analog  and  digitally  modulated  signals 
in  communication  systems.  It  can  also  identify  the  type  of  modulation.  Programs  for  Signal 
Processing  based  on  compiled  languages  such  as  FORTRAN  or  C  are  not  user  friendly,  and 
MATLAB  based  programs  have  become  the  de  facto  language  and  tools  for  signal  processing 
engineers  worldwide. 

This  thesis  describes  the  implementation  in  MATLAB  of  two  fast  methods  of  computing  the 
Spectral  Correlation  Density  (SCD)  Function  estimate,  the  FFT  Accumulation  Method  (FAM)  and  the 
Strip  Spectral  Correlation  Algorithm  (SSCA),  to  perform  Cyclic  Analysis.  Both  methods  are  based  on 
the  Fast  Fourier  Transform  (FFT)  algorithm.  The  results  are  presented  and  areas  of  possible 
enhancement  for  propeller  noise  detection  and  identification  are  discussed. 


14.  SUBJECT  TERMS 

Cyclic  Spectral  Analysis,  FFT  Accumulation  Method,  Strip  Spectral  Correlation 
Algorithm. 


15.   NUMBER  OF 
PAGES 
126 


16.   PRICE  CODE 


17.  SECURITY  CLASSIFI- 
CATION OF  REPORT 
Unclassified 


18.    SECURITY  CLASSIFI- 
CATION OF  THIS  PAGE 
Unclassified 


19.    SECURITY  CLASSIFI- 
CATION OF  ABSTRACT 
Unclassified 


20. 


LIMITATION  OF 

ABSTRACT 

UL 


NSN  7540-01-280-5500 


Standard  Form  298  (Rev.  2-89) 

Prescribed  by  ANSI  Std  239-18  298-102 


Approved  for  public  release;  distribution  is  unlimited. 

DETECTION  AND  IDENTIFICATION  OF  CYCLOSTATIONARY  SIGNALS 

Evandro  Luiz  da  Costa 

Lieutenant  Commander,  Brazilian  Navy 

B.S.,  Institute  Militarde  Engenharia,  1980 

Submitted  in  partial  fulfillment 
of  the  requirements  for  the  degree  of 

MASTER  OF  SCIENCE  IN  ELECTRICAL  ENGINEERING 

AND 
MASTER  OF  SCIENCE  IN  ENGINEERING  ACOUSTICS 

from  the 

NAVAL  POSTGRADUATE  SCHOOL 
March,  1996 


OX  LIBRARY 

'      DUATE  SCHOOL 

•Y  CA    91' 


ABSTRACT 


Propeller  noise  can  be  modeled  as  an  amplitude  modulated  (AM)  signal. 
Cyclic  Spectral  Analysis  has  been  used  successfully  to  detect  the  presence  of 
analog  and  digitally  modulated  signals  in  communication  systems.  It  can  also  identify 
the  type  of  modulation.  Programs  for  Signal  Processing  based  on  compiled 
languages  such  as  FORTRAN  or  C  are  not  user  friendly,  and  MATLAB  based 
programs  have  become  the  de  facto  language  and  tools  for  signal  processing 
engineers  worldwide. 

This  thesis  describes  the  implementation  in  MATLAB  of  two  fast  methods  of 
computing  the  Spectral  Correlation  Density  (SCD)  Function  estimate,  the  FFT 
Accumulation  Method  (FAM)  and  the  Strip  Spectral  Correlation  Algorithm  (SSCA),  to 
perform  Cyclic  Analysis.  Both  methods  are  based  on  the  Fast  Fourier  Transform 
(FFT)  algorithm.  The  results  are  presented  and  areas  of  possible  enhancement  for 
propeller  noise  detection  and  identification  are  discussed. 
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I.    INTRODUCTION 

A.  MOTIVATION 

Propeller  related  acoustic  signatures  typically  exhibit  modulation 
characteristics.  These  modulation  characteristics  originate  from  the  cavitation 
process  that  takes  place  in  the  water  due  to  the  cyclic  movement  of  the  propeller. 

The  cavitation  process  is  basically  the  collapse  of  air  and  vapor  bubbles 
due  to  variations  in  the  static  pressure.  These  variations  in  static  pressure  are  a 
consequence  of  the  passage  of  the  propeller  blades  through  the  water.  This 
movement,  cyclic  in  nature,  causes  amplitude  modulation  in  the  static  pressure, 
and  as  a  consequence  an  amplitude-modulated  (AM)  signal  can  be  detected  in  a 
receiver. 

Cyclostationary  processing  techniques  have  been  used  to  detect  and 
identify  analog  and  digital  communication  signals  very  successfully.  These 
techniques  have  the  advantage  of  using  a  more  realistic  model  for  the  signal 
than  the  stationary  model  used  in  most  of  the  more  conventional  signal 
processing  techniques. 

B.  BACKGROUND 

The  basic  elements  of  cyclic  spectral  analysis  are  the  time-variant  cyclic 
periodogram  and  the  time-variant  cyclic  correlogram.  These  two  functions  form  a 
Fourier  transform  pair.  This  fact  is  known  as  the  cyclic  Wiener  relation  or  the 
cyclic  Wiener-Khinchin  relation  [Ref.  1 :  p.  49.] 

1 


In  order  to  estimate  the  cyclic  spectrum  of  a  signal  of  interest,  both  the 
time-smoothed  cyclic  periodogram  and  the  frequency-smoothed  cyclic 
periodogram  can  be  used,  giving  rise  to  two  classes  of  computational  algorithms: 
the  time-smoothing  algorithms  and  the  frequency-smoothing  algorithms. 

Although  both  classes  of  algorithms  produce  similar  approximations  to  the 
cyclic  spectrum,  time-smoothing  algorithms  are  considered  to  be  more 
computationally  efficient  for  general  cyclic  spectral  analysis  [Ref.  2:  p. 38].  Taking 
advantage  of  the  efficiency,  two  computationally  fast  algorithms  based  on  the 
time-smoothed  cyclic  periodogram  were  developed  by  Roberts  et  al.  [Ref.  2]:  the 
FFT  Accumulation  Method  (FAM)  and  the  Strip  Spectral  Correlation  Algorithm 
(SSCA). 
C.  THESIS  GOALS 

The  purpose  of  this  thesis  is  to  implement  the  FFT  Accumulation  Method 
as  well  as  the  Strip  Spectral  Correlation  Method  in  MATLAB  [  Ref.  3].  These  two 
methods  are  already  implemented  in  C  [Ref.  4].  MATLAB  is  a  more  user-friendly 
and  widely-used  language  that  makes  simulations  and  evaluations  accessible  for 
students  and  researchers.  These  cyclic  spectral  analysis  programs  written  in 
MATLAB  can  easily  be  modified  and  are  transportable  across  operating  systems 
(i.  e.,  UNIX,  PC  systems,  MAC  systems,  VMS,  etc). 

Both  programs  are  used  to  compute  the  spectral  correlation  density  (SCD) 
function  estimate  for  a  number  of  analog  and  digital  modulated  signals,  such  as 


AM,  PAM,  ASK,  and  BPSK.  The  simulation  results  are  then  compared  to  the 
theoretical  ones.  After  that  the  attention  is  focused  on  the  AM  signals  for  which  a 
number  of  different  signal-to-noise  ratios  (SNR)  are  processed  using  both 
methods.  The  results  are  presented  and  discussed  for  two  types  of  modulation 
messages:  a  periodic  message  (a  tone)  and  a  random  message  (white  noise). 

Results  are  presented  by  showing  a  three-dimensional  plot  of  the  cyclic 
spectrum  surface,  a  contour  plot,  a  plot  of  the  power  spectral  density  (PSD) 
obtained  by  setting  the  value  of  the  cyclic  frequency  equal  to  zero,  and  some 
additional     two-dimensional      plots     for     cyclic     frequencies     of     interest. 


II.  NOISE  IN  THE  OCEAN 

Ross  [Ref.5]  defines  noise  as  unwanted  sound.  The  noise  presence 
interferes  with  the  signals  that  are  of  interest.  If  one  is  interested  in  detecting  the 
presence  of  a  particular  class  of  surface  ships,  the  sound  generated  by  a  near  by 
submarine  can  be  interpreted  as  noise.  The  reverse  situation  also  leads  to  a 
similar  statement.  Therefore,  the  concept  of  noise  has  no  absolute  definition.  It  is 
a  relative  concept,  and  its  characterization  depends  on  the  signal  of  interest  in  a 
particular  situation. 

In  this  thesis,  as  we  are  interested  in  the  detection  and  possible 
identification  of  the  noise  radiated  by  propeller,  underwater  noise  is  defined  as 
any  sound  that  interferes  with  our  ability  to  detect  and  identify  those  signals. 

As  reported  by  Ross  [Ref.  5:  p.  1],  underwater  noise  is  defined  as  sound 
in  the  water  that  limits  the  military  effectiveness  of  naval  detection  and 
identification  systems.  Noise  is  unavoidable  and  sources  that  radiate  as  much  as 
one  watt  of  acoustics  power  can  be  detected  at  relatively  long  ranges  by  modern 
passive  sonars. 

A.  TYPES  OF  UNDERWATER  NOISE 

There  are  several  different  sources  of  underwater  noise  that  are  grouped 
and  classified  in  different  ways  according  to  the  reference  used.  The  main 
sources  of  underwater  noise  can  be  divided  in  ambient  noise,  self  noise,  and 


radiated  noise.  Ambient  noise  and  self-noise  together  constitute  what  is  called 
the  background  noise  which  interferes  with  the  operation  of  a  sonar  system. 
1.  Ambient  Noise 

Ambient  noise  is  the  prevailing,  sustained  unwanted  background  of  sound 
at  some  location  in  the  ocean.  It  excludes  momentary,  occasional  sounds,  such 
as  the  noise  of  a  close-by  passage  of  a  ship  or  of  an  occasional  rain  squall.  It  is 
the  background  of  noise,  typical  of  the  location  and  depth  where  a  measuring 
hydrophone  is  located,  against  which  a  "signal,"  such  as  the  sound  of  a 
submarine  or  the  echo  from  a  target,  must  be  detected  [Ref.  6:  p.  1-1].  It 
comprises  all  noises  associated  with  the  medium  in  which  a  sonar  operates  that 
would  exist  in  the  medium  if  the  sonar  platform  or  vehicle  itself  were  not  present. 

The  spectrum  and  characteristics  of  this  kind  of  noise  are  complicated  and 
depend  upon  location,  position  of  the  receiver,  direction,  and  weather  conditions. 
In  its  most  general  form,  the  ambient  noise  spectrum  has  some  frequency  bands 
where  tonal  components  occur,  and  others  where  the  spectrum  is  continuous 
and  negatively  sloping  ("pink"  noise),  separated  by  portions  where  the  spectrum 
is  flat  ("white"  noise)  or  even  reversed  in  slope  [Ref.  6:  p.  2-1].  This  observation 
indicates  that  different  sources  of  noise  must  exist  and  be  prevalent  in  different 
regions  of  the  spectrum. 

Urick  [Ref.  6]  divides  the  overall  frequency  range  into  five  frequency 
bands:  the  ultra-low  band  (<l//z),  the  infrasonic  band  (\to20Hz),  the  low  sonic 


band  (20  to  200Hz),  the  high  sonic  band  (200to50,000//z),  and  the  ultrasonic 
{>  50kHz). 

Almost  nothing  is  known  about  the  noise  in  the  ultra-low  band,  since  just  a 
few  measurements  are  reported.  The  infrasonic  band  contains  the  strong  blade- 
rate  fundamental  frequency  of  propeller-driven  vessels,  plus  one  or  two  of  its 
harmonics.  It  is  of  great  interest  to  low  frequency  passive  sonars.  The  low  sonic 
band  is  characterized  by  the  noise  of  distant  shipping  in  areas  where  distant 
ships  are  prevalent.  In  areas  remote  from  shipping  lanes,  the  noise  in  this  band 
is  mainly  dependent  on  wind  speed  [Ref.  6].  According  to  Ross  [Ref.  5:  p.  280], 
ship-generated  noise  is  the  dominant  source  of  ambient  noise  in  the  spectral 
region  between  20  and  200  Hz  in  most  deep-water,  open  ocean  areas  and  in 
highly  traveled  seas  such  as  the  Baltic  and  Mediterranean.  The  noise  in  the  high 
sonic  band  is  very  dependent  on  the  sea  state  and  the  wind  speed.  Thermal 
noise  begins  to  dominate  the  noise  background  in  the  ultrasonic  band. 

2.  Self  Noise 

Self  noise  is  the  noise  associated  with  a  platform  and  its  sonar 
hydrophones.  It  includes  the  electronic  noise  generated  by  its  preamplifiers,  as 
observed  by  the  sonar  hydrophone  array  [Ref.  5:  p.  4].  It  can  reach  the  receiver 
by  transmission  through  the  mechanical  structure  and  by  transmission  through 
the  water  either  directly  from  the  source  or  by  reflection  from  the  sea  surface. 

Self  noise   usually  tends  to  increase  as  the  speed   of  the   platform 


increases.  At  low  frequencies  and  slow  speeds,  machinery  noise  dominates, 
whereas  at  high  frequencies  propeller  and  flow  noise  become  important. 
Actually,  as  the  speed  is  increased,  these  latter  sources  of  noise  assume  more 
importance  at  all  frequencies. 

At  very  low  speeds,  self  noise  is  usually  less  important  than  ambient 
noise.  At  higher  speeds  the  self  noise  can  become  the  limiting  factor  [Ref.  7:  p. 
413]. 

3.  Radiated  Noise 

Radiated  noise  is  that  sound  radiated  into  the  water  which  can  be  used  by 
some  receiver  to  detect  the  presence  of  the  emitter  at  a  considerable  distance. 
It  is  very  important  for  passive  sonars,  which  are  designed  to  exploit  the 
peculiarities  of  this  form  of  noise  and  to  distinguish  it  from  the  background  of  self- 
noise  and  ambient  noise  in  which  it  is  normally  observed  [Ref.  8:  p. 328]. 

A  detailed  discussion  of  the  mechanisms  involved  in  the  radiation  of 
sound  through  the  ocean  can  be  found  in  several  references  such  as  Ross  [Ref. 
5],  Kinsler  et  al.  [Ref.  7],  and  Urick  [Ref.  8]. 

B.  RADIATED  NOISE  FROM  SHIPS,  SUBMARINES,  AND 
TORPEDOES 

The  sources  of  noise  on  ships,  submarines,  and  torpedoes  can  be 
grouped    into   three   major   classes:    machinery   noise,    propeller   noise,    and 
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hydrodynamic  noise. 

Machinery  noise  comprises  that  part  of  the  total  noise  of  the  vessel 
caused  by  the  ship's  machinery.  It  originates  as  mechanical  vibration  of  the 
many  and  diverse  parts  of  a  moving  vessel.  This  vibration  is  coupled  to  the  sea 
via  the  hull  of  the  vessel. 

Machinery  vibration  can  originate  from  five  different  sources.  The  first 
source  of  machinery  noise  are  rotating  unbalanced  parts,  such  as  out-of-round 
shafts  or  motor  armatures.  The  second  one  are  repetitive  discontinuities,  such  as 
gear  teeth,  armature  slots,  and  turbine  blades.  Reciprocating  parts,  such  as  the 
explosions  in  cylinders  of  reciprocating  engines,  are  the  third  source  of 
machinery  vibration.  The  fourth  are  cavitation  and  turbulence  in  the  fluid  flow  in 
pumps,  pipes,  valves,  and  condenser  discharges.  And  mechanical  friction,  as  in 
bearings  and  journals,  is  the  fifth  source  of  machinery  noise. 

The  first  three  sources  of  machinery  vibration  produce  a  line-component 
spectrum  in  which  the  noise  is  dominated  by  tonal  components  at  the 
fundamental  frequency  and  harmonics  of  the  vibration-producing  process;  the 
other  two  give  rise  to  noise  having  a  continuous  spectrum  containing 
superimposed  line  components  from  structural  members  that  are  excited  into 
resonant  vibration.  The  machinery  noise  of  a  vessel  may  therefore  be  visualized 
as  possessing  a  low-level  continuous  spectrum  containing  strong  line 
components  that  originate  in  one  or  more  of  the  repetitive  vibration-producing 


processes  listed  above. 

Even  though  the  propeller  is  a  part  of  the  propulsion  machinery  of  a 
vessel,  the  noise  it  generates  has  both  a  different  origin  and  a  different 
frequency  spectrum  from  machinery  noise.  Machinery  noise  originates  inside  the 
vessel  and  reaches  the  water  by  various  processes  of  transmission  and 
conduction  through  the  hull.  Propeller  noise,  on  the  other  hand,  originates 
outside  the  hull  as  a  consequence  of  the  propeller  action  and  by  virtue  of  the 
vessel's  movement  through  the  water. 

The  source  of  propeller  noise  is  principally  the  cavitation  induced  by  the 
rotating  propellers  [Ref.  8:  p.  334].  Because  cavitation  noise  consists  of  a  large 
number  of  random  small  bursts  caused  by  bubble  collapse,  it  has  a  continuous 
spectrum,  covering  a  wide  frequency  range. 

Hydrodynamic  noise  originates  in  the  irregular  and  fluctuating  flow  of  fluid 
moving  by  the  vessel.  The  noise  created  by  the  turbulent  boundary  layer  is 
sometimes  called  "flow  noise."  Under  normal  circumstances,  hydrodynamic  noise 
is  likely  to  be  only  a  minor  contributor  to  radiated  noise,  and  is  apt  to  be  masked 
by  machinery  and  propeller  noises. 

C.   PROPELLER  NOISE 

Ross  [Ref.  5:  p.  253]  describes  cavitation  of  marine  propeller  as  the  most 
prevalent  source  of  underwater  sound  in  the  oceans.   Furthermore,  when  it 
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occurs,  propeller  cavitation  is  usually  the  dominant  noise  source  for  any  single 
marine  vehicle.  Submarines  and  torpedoes  often  operate  deep  enough  to  avoid 
cavitation.  Surface  ships,  on  the  other  hand,  generally  have  well-developed 
propeller  cavitation,  with  the  result  that  their  entire  radiated  spectra  from  as  low 
as  5  Hz  to  as  high  as  100  kHz  are  controlled  by  this  source  [Ref.  5:  p.  253]. 

Cavitation  is  essentially  the  rupture  of  bubbles  in  a  liquid  caused  by 
reduction  of  local  static  pressure.  It  differs  from  boiling  because  the  first  is 
caused  by  a  reduction  of  local  static  pressure  while  the  second  is  due  to  an 
increase  of  temperature.  Because  of  the  pulse  nature  of  the  individual  collapses 
and  the  random  sequence  of  occurrence,  the  resultant  spectrum  covers  a  wide 
frequency  range.  Also,  the  pulsations  of  the  aggregate  volume  of  cavitation 
radiate  strong  tonals  at  frequencies  below  70  Hz  [Ref.  5:  p.  270]. 

Of  the  various  types  of  cavitation,  blade-surface  cavitation  on  the  suction 
surface  is  the  one  that  produces  the  highest  noise  levels.  A  more  detailed 
explanation  on  the  different  types  of  cavitation,  particularly  on  the  blade-surface 
cavitation  noise  is  found  on  Ross  [Ref.  5:  pp.  253-260]. 

Radiated  spectra  of  surface  ships  are  dominated  by  propeller  cavitation 
noise  except  when  the  ships  are  operating  at  very  slow  speeds  [Ref.  5:  p.  272]. 
Some  characteristics  of  surface  ship  noise  that  confirm  the  dominance  of 
propeller  cavitation  are  strong  modulation  of  the  broadband  spectrum  at  shaft 
and  blade  frequencies  and  the  radiation  of  low-frequency  tonals  at  harmonics  of 
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these  frequencies. 

Propeller  noise  has  been  known  for  many  years  to  be  amplitude- 
modulated  and  to  contain  "propeller  beats,"  or  periodic  increases  of  amplitude, 
occurring  at  the  rotation  speed  of  the  propeller  shaft,  or  at  the  propeller  blade 
frequency,  which  is  equal  to  the  shaft  frequency  multiplied  by  the  number  of 
blades.  The  modulation  at  the  shaft  rotational  frequency  is  due  to  slight  physical 
differences  between  blades,  that  causes  one  blade  to  cavitate  more  than  the 
others.  It  is  this  shaft-rate  modulation  that  can  be  detected  by  the  human  ear  and 
which  enables  an  experienced  sonar  operator  to  determine  the  propeller 
rotational  rate  (rpm)  and  often  classify  the  target  [Ref.  5:  p.  269].  Propeller  beats 
have  long  been  used  by  listening  observers  for  target  identification  and  for 
estimating  target  speed  [Ref.  8:  p.  338]. 

Propeller  noise,  with  its  origin  in  the  flow  of  water  about  the  propeller, 
creates  tonal  components  in  addition  to  the  continuous  spectrum  of  cavitation 
noise.  Low-frequency  spectra  of  cavitating  ship  propellers  are  usually  dominated 
by  tonal  components  at  harmonics  of  the  rotational  frequency.  More  often,  at  the 
low-frequency  end  of  the  spectrum,  propeller  noise  contains  discrete  spectral 
"blade-rate"  components  occurring  at  multiples  of  the  rate  at  which  any 
irregularity  in  the  flow  pattern  into  or  about  the  propeller  is  intercepted  by  the 
propeller  blades.  The  frequency  of  the  blade-rate  series  of  line  components  is 
given  by 
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fm  =  mns,  (1) 

where  fm  is  the  frequency,  in  Hertz,  of  the  m'h  harmonic  of  the  blade-rate  series 
of  lines,  n  is  the  number  of  blades  on  the  propeller,  and  s  is  the  propeller 
rotation  speed  in  number  of  turns  per  second. 

Shaft  and  blade  modulation  frequencies  for  merchant  ships  are  now 
significantly  higher  than  they  were  during  WWII.  Forty  years  ago  most  merchant 
ships  had  three-  or  four-bladed  propellers  and  operated  at  from  60  to  100  rpm. 
Shaft  modulation  frequencies  were  generally  between  1.0  and  1.6  Hz  and  blade 
frequencies  were  from  3.5  to  6.5  Hz.  Today,  typical  merchant  propellers  have 
four,  five  or  six  blades  and  operate  at  from  75  to  135  rpm;  shaft  frequencies 
range  from  1.3  to  2.2  Hz,  and  blade  frequencies  are  typically  6  to  12  Hz  [Ref.  5: 
p.  279]. 
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III.  CYCLOSTATIONARY  PROCESSING 

The  majority  of  the  current  signal  processing  techniques  for  intercepting  or 
analyzing  manmade  signals  in  a  noise  contaminated  environment  typically 
utilize  probabilistic  models  based  on  stationary  statistics.  In  other  words,  they 
describe  the  signal  on  the  average,  and  they  have  to  restrict  themselves  to  a 
small  time  interval  in  order  for  this  approximation  to  hold.  That  limits  the  amount 
of  data  that  can  be  used  to  derive  the  features  in  the  signal  and  the  resulting 
information. 

Most  manmade  signals,  as  are  typically  encountered  in  communication, 
radar,  and  sonar  systems  have  some  parameters  that  vary  with  time.  Examples 
include  sinusoidal  carriers  in  amplitude,  phase  and  frequency  modulation 
systems,  periodic  keying  of  the  amplitude,  phase,  or  frequency  in  digital 
modulation  systems,  periodic  scanning  in  some  radar  systems,  and  amplitude 
modulation  in  propeller  noise.  This  requires  that  the  random  signal  be  modeled 
as  cyclostationary,  in  which  case  the  statistical  parameters  vary  in  time  with 
single  or  multiple  periodicity. 

Much  of  the  background  material  in  this  chapter  was  taken  from  Gardner 
[Ref.  1]. 
A.  CYCLOSTATIONARITY 

According  to  Gardner  [Ref.  1],  the  essence  of  cyclostationarity  is  the  fact 

that   sinewaves   can    be   generated   from    random   data   by   applying    certain 
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nonlinear  transformations.  As  a  consequence,  a  continuous  signal  x{t)  is 
cyclostationary  of  order  n  (in  the  wide-sense)  if  and  only  if  we  can  find  some  rih 
order  nonlinear  transformation  of  the  signal,  y{t)  =  f(x(t)) ,  that  will  generate 
finite-amplitude  additive  sine-wave  components,  which  produce  spectral  lines.  In 
the  same  sense,  a  discrete-time  signal  x[m]  is  cyclostationary  of  order  n  (in  the 

wide-sense)  if  and  only  if  we  can  find  some  rih  order  nonlinear  transformation  of 
the  signal,  y[m]  =  f[x[m]] ,  that  will  generate  finite-amplitude  additive  sine-wave 
components,  which  will  produce  spectral  lines  [Ref.  1 :  p.  2]. 

A  continuous  signal  y{t)   contains  a  finite-amplitude  additive  sine-wave 
component  with  frequency  a  ,  a  *  0,  if  the  Fourier  coefficient 

M;=(y(t)e-2^)  (2) 

is  not  zero.  In  the  same  way,  a  discrete-time  signal  y[m]  contains  a  finite- 
amplitude  additive  sine-wave  component  with  frequency  a  ,a  *  0,  if  the  Fourier 
coefficient 


,    a 
-i  ik — m 


May=[y[m)e      '•    j  (3) 

is  not  zero.  Here,  fs  denotes  the  sampling  frequency  and  i  is  the  square  root  of 
minus  one. 

The  operation  (•)  is  the  time-averaging  operation  defined  as 
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T/2 

0=lim-   \{)dt  (4) 

T — ken    /         ■» 


r/2 


for  analog  signals,  and  as 

<->=lim^r-r£(-)  (5) 

for  discrete-time  signals. 

For  second-order  cyclostationarity,  the  nonlinear  transformation  for  a 
continuous  signal  x(t)  (i.  e.,  y{t)  =  f(x(t)))  is  given  by 

yT{t)  =  x(t  +  z/2)x(t-T/2),  (6) 

while  for  a  discrete-time  signal  x[m]  (i.  e.,  y[m]  =  /[x[/w]] )  it  is  given  by 

>>Jm]  =  x[m]  x'[m-£] ,  (7) 

where  the  symbol  *  stands  for  complex  conjugation. 
B.  THE  CYCLIC  AUTOCORRELATION  FUNCTION  (ACF) 

The  Fourier  coefficient  M"  for  second-order  cyclostationarity  is  given  by 

May=  {x(t  +  r/2)x*(t-r/2)e'2nca  >.  (8) 

This  quantity  is  the  fundamental  parameter  of  second-order  periodicity  in 

continuous  time  and  is  called  the  cyclic  autocorrelation  function  (ACF),  R"(r) ,  of 

x(t). 

For  discrete-time  signals,  the  ACF  is  defined  as 
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K  M  =  (  x[n]x* [n - Z]e-ilna"  )  e'nai ,  (9) 

since  delays  other  than  sampling  increments  are  not  allowed. 

The  ACF  can  be  interpreted  as  measuring  the  amount  of  correlation 
between  different  frequency-shifted  versions  of  a  given  signal,  as  shown  in  detail 
in  Appendix  A  for  an  AM  signal. 

Therefore,  a  signal  exhibits  second-order  cyclostationarity  in  the  wide- 
sense  when  its  cyclic  autocorrelation  function,  R"(t)  for  a  continuous  time 
signal  or  R"[t]  for  a  discrete-time  signal,  is  different  from  zero  for  some  nonzero 
frequency  a  .  The  frequency  a  is  called  cycle  frequency  or  cyclic  frequency,  and 
the  discrete  set  of  cycle  frequencies  a  for  which  R°(t)*0  or  R"[£]*0  is  called 
the  cycle  spectrum  or  cyclic  spectrum. 

If  a  signal  is  cyclostationary,  the  cycle  spectrum  contains  only  harmonics 
(integer  multiples)  of  the  fundamental  cycle  frequency.  Neverthless,  if  the  signal 
has  more  than  one  fundamental  cycle  frequency,  then  the  cycle  spectrum 
contains  harmonics  of  each  of  those  frequencies.  In  this  situation  the  signal  is 
said  to  be  polycyclostationary. 
C.  THE  SPECTRAL  CORRELATION  DENSITY  FUNCTION  (SCD) 

Signals  usually  have  distinctive  features  in  the  frequency  domain  that  are 
not  easily  seen  in  the  time  domain.  Those  features  are  generally  used  for 
detecting  the  presence  of  those  signals.  For  instance,  is  very  hard  to  detect  the 
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presence  of  a  sinusoidal  signal,  when  embedded  in  noise,  by  just  looking  at  its 
time-domain  representation.  The  same  signal  can  easily  be  detected  in  the 
frequency  domain,  provided  the  integration  time  can  be  made  sufficiently  long. 

For  the  same  reason,  it  is  beneficial  to  determine  in  the  frequency  domain 
the  amount  of  correlation  between  frequency-shifted  versions  of  x{t).  The 
spectral  correlation  density  function  (SCD)  is  defined  as  the  Fourier  Transform  of 
the  cyclic  autocorrelation  function  (ACF),  to  allow  the  localization  in  the 
frequency  domain  of  the  amount  of  time-correlation  between  frequency-shifted 
versions  of  the  signal  x{t) .  The  SCD  is  given  by 

00 

S:(f)=  \Rax{r)e-2^  dr,  (10) 

-00 

for  a  continuous  signal  x(t) ,  and  by 

s:(f)=fdR:[£]e-i2afi,  (11) 

e=-ao 

for  a  discrete-time  signal  x[n] . 

To  determine  the  correlation  in  the  frequency  domain,  we  simply  pass 
both  of  the  two  frequency  translates  u{t)  =  x{t)e"nc"  and  v(t)  =  x{t)e"ra'  through 
the  same  set  of  bandpass  filters  and  then  measure  the  temporal  correlation  of 
the  filtered  signals,  as  shown  in  Figure  1  (The  signal  u{t)  represents  a  downshift 
in  frequency  by  a/2  while  v(t)  represents  an  upshift  in  frequency  by  a/2  of  the 
signal  x(t)).  In  the  limit  when  the  bandwidth  B  of  these  narrowband  components 
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approaches  zero,  we  obtain 


1 


»->o5 


S;(f)  =  )imj-(  [hi(t)®u(t)][h'B(t)®v(t)\  ), 


(12) 


where  the  symbol  8  stands  for  convolution,  and  hfB(t)  is  the  impulse  response 
of  the  bandpass  filters. 
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exp 
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Figure  1   Spectral  Correlation  Analyzer  (after  Gardner  [Ref.  1]) 


Strictly  speaking,  the  SCD  is  not  a  valid  density  function  in  the  usual 
sense,  since  it  is  not  nonnegative  and,  in  fact,  not  even  real-valued.  However, 
because  of  the  similar  properties  that  the  SCD  does  share  with  the  PSD,  the 
term  density  is  retained. 
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The  SCD  is  also  called  the  cyclic  spectral  density.  Observe  that  although 
the  argument  /  of  the  SCD  is  continuous,  as  it  always  will  be  for  a  random 
signal,  the  argument  a  is  discrete,  as  it  always  will  be  since  it  represents  the 
harmonic  frequencies  of  periodicities  underlying  the  random  time-series. 

A  detailed  example  of  the  computation  of  the  SCD  for  an  amplitude- 
modulated  (AM)  signal  is  given  in  Appendix  A. 
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IV.   ESTIMATION  OF  THE  SPECTRAL  CORRELATION  DENSITY 

FUNCTION 


Cyclic  spectral  analysis  is  used  to  detect  the  presence  of  a  signal  via  the 
spectral  correlation  density  (SCD)  function.  To  accomplish  this  goal  a  series  of 
codes  that  estimates  the  SCD  function  were  developed  in  MATLAB  language. 
Those  codes  are  implementations  of  two  FFT  based  time  smoothing  algorithms 
called  the  FFT  Accumulation  Method  (FAM)  and  the  Strip  Spectral  Correlation 
Algorithm  (SSCA).  The  majority  of  the  background  material  in  this  chapter  was 
taken  from  Roberts  et  al.  [Ref.  2]. 

As  reported  by  Roberts  et  al.  [Ref.  2],  cyclic  spectral  analysis  algorithms 
generally  fall  into  two  classes:  those  that  average  in  frequency  (frequency 
smoothing)  and  those  that  average  in  time  (time  smoothing).  Although  both 
classes  of  algorithms  produce  similar  approximations  to  the  cyclic  spectrum,  time 
smoothing  algorithms  are  considered  to  be  more  computationally  efficient  for 
general  cyclic  spectral  analysis. 

Both  methods  are  based  on  modifications  of  the  time  smoothed  cyclic 
cross  periodogram,  defined  as: 

KM/).  =lH"-/+f)  YinJ  -§))„•  (13) 

where  At  represents  the  data  time  span,  and  XT\n,f  +  —  \  and  YT\n,f-  —  \  are 

the  complex  envelopes  of  narrow-band,  bandpass  components  of  the  signals 
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x(n)  and  y{n) ,  respectively.  The  complex  envelopes  are  also  called  the  complex 
demodulates  of  the  signals.  These  complex  demodulates  are  computed  in  the 
following  way: 

r=-N'/2 

and 

YT{nj)=    I,a(r)y(n-r)e-i2*l"-')T',  (15) 

r=-N'l2 

where  a(r)\s  a  data  tapering  window  of  length  T  =  N'TS  seconds,  with  Ts  being 
the  sampling  period.  The  Fourier  transform  of  a(r)  plays  the  role  of  a  spectral 
window.  The  particular  shape  of  window,  especially  the  spectral  window,  is  of 
considerable  importance.  Windows  other  than  the  rectangle  have  a  tapering 
effect  on  the  data  they  multiply,  since  data  occurring  away  from  the  aperture 
center  are  attenuated  relative  to  the  data  at  the  aperture  center.  A  data  tapering 
window  whose  Fourier  transform  has  low  skirts  and  low  sidelobes  is  desirable  to 
reduce  cycle  leakage.  A  Hamming  window  is  used  for  the  input  bandpass  filters. 

Figure  2  shows  a  basic  implementation  of  the  discrete  time  smoothed 
cyclic  cross  periodogram,  where  the  symbol  *  stands  for  complex  conjugation. 
The  complex  demodulate  frequencies  /,  and  f2  are  related  to  the  spectrum 
frequency  /  and  the  cyclic  frequency  a  of  the  point  estimate  by  the  following 
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two  equations: 


/  = 


/i+/2 


(16) 


and 


<*  =  fi-A 


(17) 
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Figure  2  Discrete  Time  Smoothed  Cyclic  Cross  Periodogram  (after  Roberts  et  al. 
[Ref.  2]). 

A.   FFT  ACCUMULATION  METHOD  (FAM) 

In  this  method,  the  complex  demodulates  are  estimated  efficiently  by 
means  of  a  sliding    N' -point  FFT,  followed  by  a  downshift  in  frequency  to 
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baseband.  In  order  to  allow  for  an  even  more  efficient  estimation,  the  N' -point 
FFT  is  hopped  over  the  data  in  blocks  of  L  samples  (channelization).  This 
means  that  L  data  points  are  skipped  between  successive  computations  of  the 
N'  -point  FFT.  The  value  of  L  was  chosen  to  be  equal  to  N'/4  since,  according 
to  Reference  2,  p.462,  it  allows  for  a  good  compromise  between  maintaining 
computational  efficiency  and  minimizing  cycle  leakage  and  cycle  aliasing.  The 
value  of  N'  is  determined  according  to  the  desired  resolution  in  frequency  (A/) 
used  in  the  algorithm,  and  is  given  by 

N'=j^.  (18) 

N'  is  chosen  to  be  the  power  of  2  equal  to  or  larger  than  the  number 
given  by  Eq.  18  to  take  advantage  of  the  Fast  Fourier  Transform  (FFT)  algorithm 
without  making  use  of  zero-padding. 

After  the  complex  demodulates  are  computed  and  the  product  sequences 
between  each  one  of  them  and  the  complex  conjugate  of  the  others  are  formed, 
the  time  smoothing  is  accomplished  by  means  of  a  P  -point  FFT.  The  value  of  P 
is  determined  according  to  the  desired  resolution  in  cyclic  frequency  ( Aa ),  and  is 
given  by 

P  =  Tf~-  (19) 

L  Aa 

Again,   P  is  chosen  to  be  the  power  of  2  equal  to  or  larger  than  the 
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number  given  by  Eq.  19  to  take  advantage  of  the  FFT  algorithm  avoiding  zero- 
padding. 

Figure  3  illustrates  the  generation  of  the  complex  demodulates  while 
Figure  4  shows  the  implementation  of  the  FAM  method. 
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Figure  3  Computation  of  the  Complex  Demodulates 
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Figure  4  Implementation  of  the  FFT  Accumulation  Method 
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An  advantage  of  having  complex  demodulates  is  that  there  is  no  need  to 
worry  with  a  correction  factor  to  take  care  of  the  phase  shift  introduced  by 
overlap  processing.  The  last  multiplier  in  Figure  3  (i.  e.,  complex  exponential) 
provides  the  correction  to  compensate  for  the  overlap  processing  artifacts. 

The  MATLAB  codes  that  compute  and  plot  the  SCD  function  estimate 
using  the  FAM  method  are  called  AUTOFAM  and  CROSSFAM.  AUTOFAM 
computes  and  plots  the  auto  spectral  correlation  density  function  (amount  of 
correlation  between  frequency  shifted  versions  of  a  given  real  or  complex  valued 
signal)  estimate.  CROSSFAM  computes  and  plots  the  cross  spectral  correlation 
density  function  estimate  for  two  different  real  or  complex  valued  signals. 

The  inputs  required  for  these  programs  are  the  signal(s),  the  sampling 
frequency  (fs),  the  desired  frequency  resolution  (A/),  and  the  desired  resolution 
in  cyclic  frequency  (Aa).  In  the  case  of  two  signals,  the  sampling  frequencies 
must  be  the  same. 

The  programs  are  listed  in  Appendices  B  and  C. 
B.  STRIP  SPECTRAL  CORRELATION  ALGORITHM  (SSCA) 

In  the  SSCA  algorithm,  the  complex  demodulate  of  one  of  the  signals  is 
computed  in  the  same  way  it  is  done  for  the  FAM  method  (channelization).  The 
complex  demodulated  sequence  is  directly  multiplied  by  the  complex  conjugate 
of  the  other  signal.  Then,  the  resultant  signal  is  smoothed  in  time  by  means  of  an 
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N -point  FFT.  Here,  vV  is  the  total  number  of  data  points  (N  =  PL). 

Figure    5    shows    the    implementation    of   the    SSCA.    The    complex 

demodulated  component  XT{kL,fe)  is  obtained  as  shown  in  Figure  3. 

It  appears  that  the  complex  demodulate  XT\kL,fe)  is  at  a  sampling  rate 
fJL  due  to  the  introduction  of  the  decimation  factor  L ,  and  consequently  it  can 
not  be  directly  multiplied  by  x[n]  which  is  at  a  sampling  rate  fs .  However,  the 
demodulated  sequence  is  interpolated  to  match  the  sampling  rate  of  x[n]  by 
means  of  a  process  called  replication.  Replication  is  accomplished  by  holding  the 
value  of  each  complex  demodulate  sample  for  L  samples  [Ref.2]. 

The  MATLAB  codes  generated  to  compute  and  plot  the  SCD  function 
estimate  for  given  signal(s)  using  the  SSCA  method  are  called  AUTOSSCA  and 
CROSSSSCA.  AUTOSSCA  computes  and  plots  the  auto  spectral  correlation 
density  function  (amount  of  correlation  between  frequency  shifted  versions  of  a 
given  signal)  estimate.  CROSSSSCA  computes  and  plots  the  cross  spectral 
correlation  density  function  estimate  for  two  real  or  complex  valued  signals. 

The  inputs  required  for  these  routines  are  the  signal(s),  the  sampling 
frequency  (fs),  the  desired  frequency  resolution  (A/ ),  and  the  desired  resolution 
in  cyclic  frequency  {Aa).  In  the  case  of  two  signals,  the  sampling  frequencies 
must  be  the  same. 

The  codes  are  given  in  Appendices  D  and  E. 
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Figure  5  Implementation  of  the  SSCA 
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V.  EXPERIMENTAL  RESULTS 

In  this  chapter,  the  theoretical  results  of  computing  the  SCD  function  for 
some  analog  and  digital  modulated  signals  are  presented,  together  with  the 
experimental  results  obtained  as  output  from  the  programs  AUTOFAM  and 
AUTOSSCA.  An  interpretation  of  the  results  is  also  provided. 
A.  ANALOG  MODULATED  SIGNALS 

1.  Amplitude-Modulated  (AM)  Signal 

Consider  the  following  amplitude-modulated  (AM)  signal,  x[n],  given  by 

x[n]  =  a[n]p[n] ,  (20) 

where  a[n]  is  a  purely  stationary  low  pass  message  signal  with  power  spectral 
density  Sa\f) ,  and  p[n]  is  a  sinusoidal  carrier  wave  given  by 

p[n]  =  cos{27tf0n  +  <f>0).  (21) 


In  Eq.  21,  f0  =  FJFS  is  the  carrier  digital  frequency  and  cj>0  is  the  carrier 
phase.  F0  is  the  carrier  frequency  in  Hz  and  Fs  is  the  sampling  frequency  in 
Hz. 

The  SCD  function  for  this  amplitude-modulated  (AM)  sine  wave  is  given 
by 
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s:(/)= 


\s.{f  +  fo)  +  \sm(f-f.)9        a  =  0 

\sa(f)e±l2'°,  cc  =  ±2f0  (22) 

0,  otherwise. 


Details  of  the  derivation  are  given  in  Appendix  A. 

When  a[n]  is  a  tone,  some  of  the  assumptions  made  to  obtain  Eq.  22  are 
not  valid.  If  the  message  a[n]  is  a  tone  at  digital  frequency  fa  given  by 

a[n]  =  cos(27r  fan) ,  (23) 

then  it  is  necessary  to  go  through  the  same  steps  as  for  a  stationary  lowpass 
signal  with  no  spectral  lines  in  the  message's  power  spectral  density  (PSD),  as 
shown  in  Appendix  A. 

Let  us  assume  that   x[n]    is  an  amplitude-modulated  single-sideband 
(AMSSB)  signal,  obtained  by  transmitting  only  the  lower-side  frequency,  given  by 

x[n]=  -[cos(2tt  fan)  ■  cos(2n  f0n  +  <f>0)  +  sm{2n  fan)  ■  sin(2x  f0n  +  #,)] 

=  |cos[2;r(/0 -/>  +  &].  (24) 

Since  the  cyclic  autocorrelation  function  (ACF)  is  given  by 

R?U]  =  (x[n]  x'[n-l]  e^2*"")  eixai ,  (25) 

replacing  x[n]  and  x[n-£]  by  Eq.  24  leads  to, 

RM=~MHfo  -fa)t]e'*a((e-'2™)  + 
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]_e,n\a-2{f0-fa)\t  I    i2^[«-2(/0 -/„)]»  \     i2*0   +  }_  g»r[«+2(/0-/e)W    -/2*[«+2(/„ -/.)]»  \     -,74,  _       (26) 


This  can  be  rewritten  as 


CM= 


jcos[2^(/0-/0)4  a  =  0 

a  =  ±2(f0-fa). 


-e±a*° 
8 


(27) 


The  SCD  is  the  Fourier  transform  of  the  cyclic  autocorrelation  function. 
Therefore,  the  SCD  is  given  by  the  Fourier  transform  of  Eq.  27  leading  to 


s:(/)= 


\[s{f-fo+fahS{f  +  fo-fa)\  «  =  0 


fay** 


a  =  ±2(f0-fa). 


(28) 


Now,   let  us  suppose  that   x[n]      is  an  amplitude-modulated   double- 
sideband  supressed  carrier  (AMDSB-SC)  signal  given  by 

x[n\  =  cosllir  fan)  ■  cos[27rf0n  +  <f>0).  (29) 

This  can  be  written  using  a  trigonometric  identity  as 

x[n]  =  « {cos[2^(/0  +  fa)n  +  <j>0 ]  +  cos[2;r(/0  -  fa)n  +  ^ ]} .  (30) 

So,  replacing  x[«]  and  x[«-^]  by  Eq.  30  into  Eq.  25  leads  to, 


*;[*]=  i[cos2*(/0  +  /a)*  +  cos2*{/.  -/.)/]  e^2™) 


+ 


J_     ,>[a-2(/0+/a)]W    -i2«(«-2(/.+/,)]"\     «A     ,J_     "r[a+2(/„-/jW    -,2*(a-2/>\ 

16  \  /  16  \  I 

]_     m\a-2{f0-fa))t  I    -i2x{a-2f„)n\     ,2^     .^_Ma+2(f0+faj\el-i2n\a+2{f0+f0)\n\l2to 


16 


e  Yo  + — e 
16 


e     r°  + 
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J_fcr[«+2(/„ -/„)]*/    -i2*{a+2/a)n\     -12*.        J_     m\a-l{j 0- f 0)\t  /  ^-i2„{a+2/a)n 
16  \  /  16 


J_eM«+2(/0+/JW    -««(a+2/.)i.\       J_    /4«-2(/0 +/„)]«/    -/2*(«-2/0)«\     ,2& 

16  \  /     16  \  / 

J_     «[a-2(/„-/.)]<  /    -i2*[a-2(/e -/„)]«  \     ,2^        J_     .*[  a +2  (/„+/„)  W    ->2^(«+2/„)«  \     -,2& 

16  \  /  16  \  I 

This  can  be  rewritten  as 


*;m= 


^[cos2;r(/0  +  /.)/  +  cos2;r(/0  -  /a)4  a  =  0 

a  =  ±2/. 
«  =  ±2(/0-/J 

«  =  ±2/. 
«  =  ±2(/0+/0). 


16  16  8  ° 

— e±<2*  , 
16 

8  7a 

1 


(32) 


16 


Therefore,  the  SCD  is  given  by  the  Fourier  transform  of  Eq.  32  leading  to 


s°Af)  = 


16 


a  =  ±2(/0-/J 
«  =  ±2/0 


(33) 


j^M*"*. 


16 


K/-/*M/+/J]^ 


£<*"*. 


a  =  ±2(/0+/J. 


For  the  case  where  x[n]     is  an  amplitude-modulated  double-sideband 


transmitted  carrier  (AMDSB-TC)  signal  given  by 

x[n]  =  [l  +  cos(2;r./>)]  •  cos(2a  f0n  +  <f>0). 


(34) 
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This  can  be  written  using  a  trigonometric  identity  as 


x[n]  =  cos(2tt  f0n)  +  - {cos[2;r(/0  +  fa )n  +  fa ]  +  cos[2^(/0  -fa)n  +  <p0 ]} .  (35) 

So,  replacing  x[n]  and  x[«-^]  by  Eq.  35  into  Eq.  25  leads  to, 
RZ[t]=\\co!(2xfj)  +  Ucos24f0  +fa)i  +  cos24f0  -/.)/]}  e"  («-«--)  + 


8  8  8  8  J\  / 


+ 


II  /«[«+2(/.-/.)]<  + 1  M«-2(/.+/.)]<  +Ie*K«-2/.)<  .1  g*K«+2/.)l/ g-iM«-/>\j 
[8  8  8  8  J\ 

J_    *r(a-2{/.-/.)]/       J_    **(«+2(/#+/.)]*l/    -i2«r(a+2/„)«  \ 

,16  16  J\  / 


J_    i«[«+2(/.-/.)J<       J_    l«(a-2(/.+/.)]/\//>-/2«(a-2/.)ii 
16  16 


J_e'>[a+2(/0-/a)]^-/2^/e-/24a+2(/<,-/.)]«') 


_LgM«-2(/.-/.)]<gnA/g-'2«|«-2(/.-/.)]»,; 


fle««(«+2(/.-/.)l<  +  Ie'*(«+2/0)*  1    -nA  /  -n«{«+(2/.-/.)]«\  + 


IeM«-2(/.-/J]'  +Ie*K«-V.)Ugl2A  Ie-i2f\a-{2f0-fa)\n 


}_    i*(a+2ft)t  J_'-«|«+2(/0 -/„)]/       J_it[«+2(/,, +/„)]<  L^2A  l-iln(a+lJB)n 

,4  16                          16                      J          \ 

flgM«-2/.)'  ,_L-**{«-2(/.-/.)l<  +_Le"r[a-2U+'*)]'  L''2^/e-'"2T(«-2/.)« 

U  16                          16                      J         \ 


35 


+ 


J.     b^a+2[f.*f.)]l       }_     i*{a+2J0)t  \     -/2*  /    -i2x[a+{2fB  +/.)]« 


8 


8 


+ 


}_     i*[a-2{f0+/a)\l       }_     ix{a-2fe)l}     12*  /p-i2«[a^2/.+/.)]» 

8  8  J 


J_    /«-[«+2(/0+/a)]^    -/2A,  /_ -n«{«+2(/-+/.)]« 


16 


e     ™  \  e 


+ 


J_    i*{a-2{fc+fa)]l    njj    -i2x[a-2(f0+fa)]n 
16  \ 


(36) 


This  can  be  rewritten  as 


r:[£]  = 


\co{K(2f0  -  fa)t]  +  ^cos[^(2/0  +  /a)4 
±cw(2*/.j), 


,±'2  A, 


16 


1 


-co^7ifj)e±l 

4     8      v        a  ' 
1 


,±'2^ 


cos^/^)^, 


—e±i2*° 
16 


«  =  0 

«  =  ±2/. 
a  =  ±2(/o-/a)(37) 

a  =  ±(2/.-/.) 

a  =  ±2/0 

«  =  ±(2/0+/J 
a  =  ±2(/c+/a). 


Therefore,  the  SCD  is  given  by 


s;(fh 


2)     8 


£<c/+/.)+i*/)+£*/-/.) 


«"*, 


16 


</)« 


«-*/.    (38) 

a  =  ±2(/„-/,) 
a  =  ±(2/„-/„) 

a  =  M2f„-f.) 

a  =  ±2(f0+f.). 
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The  results  obtained  from  AUTOFAM  and  AUTOSSCA  are  presented  as  a 
sequence  of  three  plots,  for  each  of  the  signals  utilized.  This  presentation  is 
maintained  throughout  (i.  e.,  Figure  6  -  Figure  41).  The  first  plot  is  a  surface  plot 
that  shows  the  magnitude  of  the  SCD  function  estimate  in  each  region  of  the 
bifrequency  plane  with  coordinates  /  and  a.  The  second  plot  is  a  two- 
dimensional  contour  plot.  It  gives  a  better  view  of  the  position  of  the  features  in 
the  bifrequency  plane.  The  third  plot  is  a  set  of  two-dimensional  slices  of  the 
SCD  function  estimate  for  fixed  values  of  the  cyclic  frequency  a  . 

Figures  6-11  show  the  results  obtained  by  using  the  programs 
AUTOFAM  and  AUTOSSCA  on  an  amplitude-modulated  single-side  band 
(AMSSB)  signal.  For  fa  =  5\2Hz  and  f0  =  204SHz,  Equation  28  leads  to  the 
following  result: 


s;(/)= 


-[S(f- 1536)  +  S(f  + 1536)1  a  =  0 

8  i  (39) 

S(f)e±l2'°,  a  =  ±3072Hz. 


So,  according  to  Eq.  39,  ones  expect  to  obtain  peaks  at  /  =  ±\536Hz  for 
a  =  0,  and  at  /  =  0  for  a  =  ±3072Hz.  The  results  in  Figures  6-11  are  in 
agreement  with  the  theoretical  results  in  Eq.  39  and  show  that  AUTOFAM  gives 
a  clearer  picture  than  AUTOSSCA.  This  is  even  true  when  AUTOSSCA  uses  a 
finer  resolution  in  cyclic  frequency  (i.  e.,  Aa  =  64  Hz  for  AUTOFAM  versus 
Aa  =  32  Hz  for  AUTOSSCA). 

The    results    obtained    for   an    amplitude-modulated    double-side    band 
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suppressed-carrier    (AMDSB-SC)     are    presented     in     Figures     12-17.     For 
fa=5\2Hz  and  f0  =  2048Hz ,  Equation  33  leads  to  the  following  result: 


sM  = 


— \s(f  -2560) +S(f  - 1536)  +  S(f  + 1536)+S(f +2560)L  a  =  0 

16l 

[S(f  -  2048) + S(f  +  2048)], 

1 


1 
l6l 


16 


to 


~[S(/-512)+S(/+512)], 


16 


16 


4A 


a  =  ±\024Hz 
a  =  ±3072/£ 
a  =  ±4096Hz 
a  =  ±5\20Hz. 


(40) 


According  to  Eq.  40,  ones  expect  to  have  four  peaks  at  /  =  ±\536Hz  and 
/  =  ±2560Hz ,  for  a  =  0;  two  peaks  at  /  =  ±204SHz ,  for  a  =  ±\024Hz ;  one  peak 
at  /  =  0 ,  for  a  =  ±3072Hz ;  two  peaks  at  /  =  ±5\2Hz ,  for  a  =  ±4096#z ;  and  a 
peak  at  /  =  0 ,  for  a  =  ±5\20Hz .  Figures  12-17  confirm  the  theoretical  results  as 
given  by  Eq.  40. 

The  results  for  an  amplitude-modulated  double-side  band  transmitted- 
carrier  (AMDSB-TC)  are  presented  in  Figures  18-23.  For  fa=5\2Hz  and 
f0  =  2048  Hz ,  Equation  38  leads  to  the  following  result: 


*?(/)= 


—  S{/  +  2560)  +  -S(/  +  2048)  +  —  S(f  + 1536)  +  —  <5{/- 1536)  +  -  S(f  -  2048)  +  —  S(f  -  2560),  a  =  0 

-  S(f  +  2304)  +  -  S(f  + 1 792)  +  -  S(f  - 1 792)  +  -  S{f  -  2304), 

o  o  o  o 

!<5{/  +  2048)  +  -i5(/-2048). 


\b 


w*-. 


S{/  +  256)  +  -S(f  -256) 
I<$(/  +  256)  +  I<5(/-256) 


bw 


a  =  ±512* 

a  =  ±1024//r 
a  =  ±3072/£ 
a  =  ±3584/fc 
a  =  ±4096//: 
a  =  ±4608//r 
a  =  ±5120//r. 


(41) 


According  to  Eq.  41,  ones  expect  to  have  two  big  peaks  at  /  =  ±204$Hz 
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and  four  smaller  peaks  at  /  =  ±\536Hz  and  /  =  ±2560 Hz ,  for  a  =  0;  four  peaks 
at  /  =  ±\792Hz  and  /  =  ±2304 Hz  ,  for  a  =  ±5\2Hz\  two  peaks  at  /  =  ±2048 Hz, 
for  a  =  ±1024//r ;  one  peak  at  /  =  0 ,  for  a  =  ±3072//- ;  two  peaks  at  /  =  ±256Hz , 
for  a  =  ±3584//r ;  a  big  peak  at  /  =  0  and  two  smaller  peaks  at  /  =  ±512//z  ,  for 
a  =  ±4096Hz  ;  two  peaks  at  /  =  ±256Hz ,  for  a  =  ±4608//z  ;  and  a  peak  at  /  =  0 , 
for  a  =  ±5120//r  .  Figures  18-23  verify  the  theoretical  results  of  Eq.  41. 

To  have  a  reliable  estimation  of  the  spectral  correlation  density  function  is 
necessary  that  the  product  A/A/  »1  [Ref.  1].  This  condition  requires  that 
A/»  Aa .  In  some  of  the  results  obtained  (i.  e.,  as  in  Figure  21),  this 
requirement  was  not  met.  Computational  limitations  did  not  allow  for  a  better 
resolution  in  the  plots,  since  this  translates  to  generating  a  large  amount  of  data, 
making  impossible  to  obtain  a  plot  on  the  printer  and/or  on  the  screen  with  the 
available  PC  and  workstation  resources. 
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Figure  6  Surface  plot  of  the  SCD  estimate  magnitude  for  an  Amplitude 
Modulated  (AM)  Single-Side  Band  Signal,  using  AUTOFAM,  with  the  following 
parameters:  &f  =  256Hz ,  Act  =  64 Hz ,  fc  =  2048//z ,  /  =5\2Hz ,  and  fs  =  8192/fe , 
where  fc  is  the  carrier  frequency,  /,  is  the  tonal  frequency,  and  fs  is  the 
sampling  frequency. 
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Figure  7  Contour  plot  of  the  SCD  estimate  magnitude  for  an  Amplitude 
Modulated  (AM)  Single-Side  Band  Signal,  using  AUTOFAM,  with  the  following 
parameters:  Af  =  256Hz,  &a  =  64Hz,  fc=204SHz,  ft  =512/fe,and  fs  =  S\92Hz, 
where  fc  is  the  carrier  frequency,  /,  is  the  tonal  frequency,  and  fs  is  the 
sampling  frequency. 
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Figure  8  Plots  of  the  SCD  estimate  magnitude  for  an  Amplitude  Modulated  (AM) 
Single-Side  Band  Signal,  using  AUTOFAM,  for  a  =  0  and  a  =  3072Hz, 
respectively,  with  the  following  parameters:  Af  =  256Hz,  Aa  =  64Hz, 
fc  =  2048Hz ,  /,  =  5\2Hz ,  and  fs  =  %\92Hz ,  where  fc  is  the  carrier  frequency,  /, 
is  the  tonal  frequency,  and  f  is  the  sampling  frequency. 
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Figure  9  Surface  plot  of  the  SCD  estimate  magnitude  for  an  Amplitude 
Modulated  (AM)  Single-Side  Band  Signal,  using  AUTOSSCA,  with  the  following 
parameters:  Af  =  256Hz,  Aa  =  32Hz,  /c=2048tfz,  f,=512Hz,  and  fs  =  S\92Hz, 
where  fc  is  the  carrier  frequency,  /,  is  the  tonal  frequency,  and  fs  is  the 
sampling  frequency. 
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Figure  10  Contour  plot  of  the  SCD  estimate  magnitude  for  an  Amplitude 
Modulated  (AM)  Single-Side  Band  Signal,  using  AUTOSSCA,  with  the  following 
parameters:  Af  =  256Hz,  Aa  =  32Hz,  /c=2048//z,  ft=5\2Hz,  and  fs  =  &\92Hz, 
where  fc  is  the  carrier  frequency,  /,  is  the  tonal  frequency,  and  fs  is  the 
sampling  frequency 
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Figure  1 1  Plots  of  the  SCD  estimate  magnitude  for  an  Amplitude  Modulated 
(AM)  Single-Side  Band  Transmitted  Carrier  Signal,  using  AUTOSSCA,  for  a  =  0 
and  a  =  3072Hz,  respectively,  with  the  following  parameters:  Af  =  256Hz, 
Aa  =  32Hz,  fc  =  2048Hz,  ft=5l2Hz,  and  fs  =  S\92Hz,  where  fc  is  the  carrier 
frequency,  ft  is  the  tonal  frequency,  and  fs  is  the  sampling  frequency. 
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Figure  12  Surface  plot  of  the  SCD  estimate  magnitude  for  an  Amplitude 
Modulated  (AM)  Double-Side  Band  Supressed  Carrier  Signal,  using  AUTOFAM, 
with  the  following  parameters:  A/  =  256/fe ,  Aa  =  32Hz ,  fc  =  2048Hz ,  ft  =5\2Hz , 
and  fs  =  8\92Hz,  where  fc  is  the  carrier  frequency,  ft  is  the  tonal  frequency, 
and  /  is  the  sampling  frequency. 
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Figure  13  Contour  plot  of  the  SCD  estimate  magnitude  for  an  Amplitude 
Modulated  (AM)  Double-Sided  Band  Supressed  Carrier  Signal,  using 
AUTOFAM,  with  the  following  parameters:  A/  =  256Hz ,  Aa  =  32Hz ,  fc  =  204SHz , 
ft  =512Hz,  and  fs  =  8\92Hz,  where  fc  is  the  carrier  frequency,  /,  is  the  tonal 
frequency,  and  fs  is  the  sampling  frequency. 
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Figure  14  Plots  of  the  SCD  estimate  magnitude  for  an  Amplitude  Modulated 
(AM)  Double-Side  Band  Supressed  Carrier  Signal,  using  AUTOFAM,  for  a  =  0 , 
a  =  1024Hz,  a  =  3072Hz,  a  =  4096Hz,  and  a  =  5\20Hz,  respectively,  with  the 
following  parameters:  A/  =  256Hz ,  Act  =  32Hz ,  fc  =  204SHz ,  f,=5\2Hz,  and 
/,  =  81927/z ,  where  fe  is  the  carrier  frequency,  /,  is  the  tonal  frequency,  and  fs 

is  the  sampling  frequency. 

48 


alpha(Hz) 


Figure  15  Surface  plot  of  the  SCD  estimate  magnitude  for  an  Amplitude 
Modulated  (AM)  Double-Side  Band  Supressed  Carrier  Signal,  using 
AUTOSSCA,  with  the  following  parameters:  Af  =  256Hz,  Aa  =  32Hz, 
fc  =  2048Hz ,  f,  =  5\2Hz ,  and  /,  =  S\92Hz ,  where  fc  is  the  carrier  frequency,  ft 
is  the  tonal  frequency,  and  /  is  the  sampling  frequency. 
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Figure  16  Contour  plot  of  the  SCD  estimate  magnitude  for  an  Amplitude 
Modulated  (AM)  Double-Sided  Band  Supressed  Carrier  Signal,  using 
AUTOSSCA,  with  the  following  parameters:  Af  =  256Hz,  Aa  =  32Hz, 
fc  =  2048#z ,  ft  =  5\2Hz ,  and  fs  =  S\92Hz ,  where  fc  is  the  carrier  frequency,  ft 
is  the  tonal  frequency,  and  /  is  the  sampling  frequency. 
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Figure  17  Plots  of  the  SCD  estimate  magnitude  for  an  Amplitude  Modulated 
(AM)  Double-Side  Band  Supressed  Carrier  Signal,  using  AUTOSSCA,  for  a  =  0 , 
a  =  \024Hz,  a  =  3072Hz,  a  =  4096Hz,  and  a  =  5\20Hz,  respectively,  with  the 
following  parameters:  Af  =  256Hz,  Aa  =  32Hz,  fc=204SHz,  ft=5\2Hz,  and 
fs  =  Sl92Hz ,  where  fc  is  the  carrier  frequency,  ft  is  the  tonal  frequency,  and  fs 

is  the  sampling  frequency. 
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Figure  18  Surface  plot  of  the  SCD  estimate  magnitude  for  an  Amplitude 
Modulated  (AM)  Double-Side  Band  Transmitted  Carrier  Signal,  using  AUTOFAM, 
with  the  following  parameters:  A/  =  128/fe ,  Act  =  32Hz ,  fe  =  2048//z ,  /,  =  512Hz , 
and  fs  =  S\92Hz,  where  fc  is  the  carrier  frequency,  ft  is  the  tonal  frequency, 
and  f  is  the  sampling  frequency. 
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Figure  19  Contour  plot  of  the  SCD  estimate  magnitude  for  an  Amplitude 
Modulated  (AM)  Double-Sided  Band  Transmitted  Carrier  Signal,  using 
AUTOFAM,  with  the  following  parameters:  A/  =  128/fe ,  Aa  =  32Hz ,  fc  =  2048/fz , 
ft  =5\2Hz,  and  fs  =  8\92Hz,  where  fc  is  the  carrier  frequency,  ft  is  the  tonal 
frequency,  and  /  is  the  sampling  frequency. 
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Figure  20  Plots  of  the  SCD  estimate  magnitude  for  an  Amplitude  Modulated 
(AM)  Double-Side  Band  Transmitted  Carrier  Signal,  using  AUTOFAM,  for  a  =  0, 
a  =  \024Hz,  a  =  3072Hz,  a  =  4096Hz,  and  a  =  5\20Hz,  respectively,  with  the 
following  parameters:  Af  =  l2SHz,  Aa  =  32Hz,  fc=2048Hz,  ft=5\2Hz,  and 
fs  =  S\92Hz ,  where  fc  is  the  carrier  frequency,  /,  is  the  tonal  frequency,  and  fs 

is  the  sampling  frequency. 
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Figure  21  Surface  plot  of  the  SCD  estimate  magnitude  for  an  Amplitude 
Modulated  (AM)  Double-Side  Band  Transmitted  Carrier  Signal,  using 
AUTOSSCA,  with  the  following  parameters:  Af  =  \2SHz,  Aa  =  128Hz, 
fc=204SHz,  ft  =512Hz,  and  fs  =  S\92Hz,  where  fc  is  the  carrier  frequency,  /, 
is  the  tonal  frequency,  and  f  is  the  sampling  frequency. 
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Figure  22  Contour  plot  of  the  SCD  estimate  magnitude  for  an  Amplitude 
Modulated  (AM)  Double-Sided  Band  Transmitted  Carrier  Signal,  using 
AUTOSSCA,  with  the  following  parameters:  Af  =  \2SHz,  Aa  =  32Hz, 
fc  =  2048//z ,  /,  =5\2Hz ,  and  fs  =  8192//z ,  where  fc  is  the  carrier  frequency,  /, 
is  the  tonal  frequency,  and  fs  is  the  sampling  frequency.     . 
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Figure  23  Plots  of  the  SCD  estimate  magnitude  for  an  Amplitude  Modulated 
(AM)  Double-Side  Band  Transmitted  Carrier  Signal,  using  AUTOSSCA,  for 
a  =  0,  a  =  \024Hz,  a  =  3072Hz,  a  =  4096Hz,  and  a  =  5120Hz,  respectively,  with 
the  following  parameters:  Af  =  32Hz,  Aa  =  32Hz,  fc  =  204SHz,  /,  =512ife,and 
/,  =  S\92Hz ,  where  fe  is  the  carrier  frequency,  /,  is  the  tonal  frequency,  and  /, 

is  the  sampling  frequency. 
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2.  Pulse-Amplitude  Modulated  (PAM)  Signal 

In  the  previous  section,  consider  the  following  form  for  the  carrier  wave, 

00 

p[n]=  ^S(n-mN0)  (42) 

OT=-oo 

where  N0  is  the  pulse  period.  If  the  product  x[n]  =  a[n]p[n]  is  filtered  using  a 
pulse  form  with  impulse-response  q[n],  then  the  result  is  the  pulse-amplitude 
modulation(PAM)  signal 

00 

y[n]=  X  a[mNo  ]  l[n  ~  mN o  ]  •  (43) 

m=-oo 

The  cyclic  spectra  for  the  PAM  signal,  when  a[n]  is  stationary,  is  given  by 


s;{/)= 


T. 


a  = 


k 
T0 


where  Q{f)  is  the  Fourier  transform  of  the  pulse  form  q[n],  and  Sa(f)  is  the 
power  spectral  density  of  the  signal  a[n] . 

Figures  24-29  show  surface  and  contour  plots  of  the  SCD  estimate 
magnitude  for  a  PAM  signal  which  is  modulated  by  a  random  sequence  of  zeros 
and  ones. 

From  Eq.  44,  ones  expect  to  obtain  something  different  from  zero  just  at 
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cyclic  frequencies  that  are  multiple  of  the  bit  rate.  Figures  24-29  verify  the  result 
obtained  in  Eq.  44.  The  poor  resolution  is  more  obvious  in  these  examples, 
making  it  difficult  to  obtain  the  appropriate  SCD  representation  for  the  signal. 
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Figure  24  Surface  plot  of  the  SCD  estimate  magnitude  for  a  Pulse-Amplitude 
Modulated  (PAM)  Signal,  using  AUTOFAM,  with  the  following  parameters: 
Af  =  5\2Hz,  Aa  =  \6Hz,  rb=\024Hz,  and  fs=S\92Hz,  where  rb  is  the  bit  rate, 
and  f  is  the  sampling  frequency. 
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Figure  25  Contour  plot  of  the  SCD  estimate  magnitude  for  a  Pulse-Amplitude 
Modulated  (PAM)  Signal,  using  AUTOFAM,  with  the  following  parameters: 
Af  =  5l2Hz,  Aa  =  \6Hz,  rb  =  \024Hz,  and  fs=S\92Hz,  where  rb  is  the  bit  rate, 
and  fs  is  the  sampling  frequency. 
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Figure  26  Plots  of  the  SCD  estimate  magnitude  for  a  Pulse-Amplitude 
Modulated  (PAM)  Signal,  using  AUTOFAM,  for  a  =  0,  a  =  1024//-,  a  =  204SHz, 
a  =  3072Hz,  and  a  =  4096Hz,  respectively,  with  the  following  parameters: 
A/  =  512//z ,  Aa  =  \6Hz  ,  rb=  \024Hz ,  and  /,  =S\92Hz ,  where  rb  is  the  bit  rate, 
and  fs  is  the  sampling  frequency. 
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Figure  27  Surface  plot  of  the  SCD  estimate  magnitude  for  a  Pulse-Amplitude 
Modulated  (PAM)  Signal,  using  AUTOSSCA,  with  the  following  parameters: 
Af  =  5l2Hz,  Aa  =  \6Hz,  rb=\024Hz,  and  fs  =  8\92Hz,  where  rb  is  the  bit  rate, 
and  fs  is  the  sampling  frequency. 
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Figure  28  Contour  plot  of  the  SCD  estimate  magnitude  for  a  Pulse-Amplitude 
Modulated  (PAM)  Signal,  using  AUTOSSCA,  with  the  following  parameters: 
Af  =  5\2Hz,  Aa  =  \6Hz,  rb=\024Hz,  and  fs  =  S\92Hz,  where  rh  is  the  bit  rate, 
and  fs  is  the  sampling  frequency. 
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Figure  29  Plots  of  the  SCD  estimate  magnitude  for  a  Pulse-Amplitude 
Modulated  (PAM)  Signal,  using  AUTOSSCA,  for  a  =  0 ,  a  =  1024Hz ,  a  =  2048Hz , 
a  =  3072Hz,  and  a  =  4096Hz,  respectively,  with  the  following  parameters: 
Af  =  5\2Hz,  Aa  =  \6Hz,  rb=\024Hz,  and  fs  =  $192Hz,  where  rb  is  the  bit  rate, 
and  f  is  the  sampling  frequency. 
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B.    DIGITAL  MODULATED  SIGNALS 

1.    Amplitude  Shift  Keying  (ASK)  Signal 

An  ASK  signal  is  simply  an  AM  signal 


( 


fo 


\ 


x[n]  =  a[n]  cos  27r-f~n  +  0o 

v        fs 


in  which  the  amplitude  message  a[n]  is  a  M-ary  PAM  signal 


(45) 


a[n]=  Y,ak  q{n-kN0-n0). 

k=-x 

The  spectral  correlation  density  function  for  this  ASK  signal  is  given  by 


(46) 


s^=M 


a 


Q\f  +  fo  +  2J 


Q'{f+f0-j)s;(f+f0) 


+  Q{f-f0+^]Q,{f-f0-fjs:(f-f0) 


-ilnan. 


+ fl(/ + § + /.)  Q'{f  -  §  -  /.)  srf-  (/)  e-'N«^>.^.i 

+ e(/ + f  -  /.)  e'(/  -  § + /.)  sr!/-  (/)  e-M"-^-2*i } .      (47) 


For  a  full-duty  cycle  rectangle  pulse,  q[n]  is  given  by 


q[n]  = 


X     \n\<Njl 
0,    H>iV./2 


Its  Fourier  transform  is  given  by 


(48) 
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sin(;r>2Vj 

Q{f)  =  ~r •  (49) 

Details  on  this  result  can  be  found  in  Reference  9. 

Figures  30-41  present  the  outputs  from  AUTOFAM  and  AUTOSSCA  for 
some  ASK  signals.  Figures  30-35  present  the  results  for  an  ASK  signal  with  bit 
rate  rb  =  2048Hz ,  and  Figures  36-41  present  the  results  for  the  same  signal  with 
rb=\024Hz.  As  f0  =  2048Hz ,  from  Eq.  47,  ones  expect  to  obtain  a 
representation  which  is  a  combination  of  four  PAM  signals  centered  at 
f  =  ±fo=  ±204SHz ,  for  a  =  0,  and  at  /  =  0 ,  for  a  =  ±2f0  =  ±4096Hz .  The  results 
obtained  in  Figures  30-41  confirm  that  expectation. 

2.    Binary-Phase  Shift  Keying  (BPSK)  Signal 

A  phase-shift  keying  (PSK)  signal  is  just  a  phase-modulated  (PM)  carrier 

(  f0      \ 

x[n]  =  cos  2n— r«  +  ^[«]  (50) 

in  which  the  phase  <p[n]  is  a  M-ary  PAM  signal 

00 

4»]=  Y,a\mNo]An-mNo]-  (51) 

Hence,  the  cyclic  spectra  for  the  BPSK  signal,  x[n],  is  given  by  Eq.  47, 
already  computed  for  a  M-ary  ASK  signal.  Of  course,  the  results  are  the  same 
as  for  the  M-ary  ASK  signal  and  are  given  by  Figures  30-41 . 
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Figure  30  Surface  plot  of  the  SCD  estimate  magnitude  for  an  Amplitude  Shift 
Keying  (ASK)  or  a  Binary  Phase  Shift  Keying  (BPSK)  Signal,  using  AUTOFAM, 
with  the  following  parameters:  A/  =  512/fe ,  Act  =  \6Hz ,  fc  =  2048ife , 
rh=204$Hz,  and  fs  =  8\92Hz,  where  fc  is  the  carrier  frequency,  rb  is  the  bit 
rate,  and  fs  is  the  sampling  frequency. 
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Figure  31  Contour  plot  of  the  SCD  estimate  magnitude  for  an  Amplitude  Shift 
Keying  (ASK)  or  a  Binary  Phase  Shift  Keying  (BPSK)  Signal,  using  AUTOFAM, 
with  the  following  parameters:  Af  =  5l2Hz,  Aa  =  16Hz,  fc=2048Hz, 
rb=2048Hz,  and  fs  =  8\92Hz,  where  fc  is  the  carrier  frequency,  rb  is  the  bit 
rate,  and  /  is  the  sampling  frequency. 


8000 


69 


ft 


-5000  -4000  -3000  -2000  -1000 


0 
f(Hz) 


1000      2000      3000      4000 


Figure  32  Plots  of  the  SCD  estimate  magnitude  for  an  Amplitude  Shift  Keying 
(ASK)  or  a  Binary  Phase  Shift  Keying  (BPSK)  Signal,  using  AUTOFAM,  for 
a  =  0 ,  a  =  204SHz ,  a  =  4096Hz ,  and  a  =  6144 Hz  ,  respectively,  with  the  following 
parameters:  Af  =  5\2Hz,  Aa  =  \6Hz,  fc=204SHz,  rb=2048Hz,  and 
fs  =  $\92Hz,  where  fc  is  the  carrier  frequency,  rb  is  the  bit  rate,  and  fs  is  the 

sampling  frequency. 

70 


alpha(Hz) 


Figure  33  Surface  plot  of  the  SCD  estimate  magnitude  for  an  Amplitude  Shift 
Keying  (ASK)  or  a  Binary  Phase  Shift  Keying  (BPSK)  Signal,  using  AUTOSSCA, 
with  the  following  parameters:  A/  =  1024Hz,  Aa  =  8Hz,  fc  =  204SHz, 
rb=204SHz,  and  fs  =  S\92Hz,  where  fc  is  the  carrier  frequency,  rb  is  the  bit 
rate,  and  fs  is  the  sampling  frequency. 
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Figure  34  Contour  plot  of  the  SCD  estimate  magnitude  for  an  Amplitude  Shift 
Keying  (ASK)  or  a  Binary  Phase  Shift  Keying  (BPSK)  Signal,  using  AUTOSSCA, 
with  the  following  parameters:  Af  =  \024Hz,  Aa  =  SHz,  fc=204SHz, 
r6=2048/fe,  and  fs=8\92Hz,  where  fe  is  the  carrier  frequency,  rb  is  the  bit 
rate,  and  /  is  the  sampling  frequency. 
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Figure  35  Plots  of  the  SCD  estimate  magnitude  for  an  Amplitude  Shift  Keying 
(ASK)  or  a  Binary  Phase  Shift  Keying  (BPSK)  Signal,  using  AUTOSSCA,  for 
a  =  0 ,  a  =  2048/fe ,  a  =  4096Hz ,  and  a  =  6U4Hz ,  respectively,  with  the  following 
parameters:  A/  =  1024#z,  Aa  =  SHz,  fc=204SHz,  rb=  2048Hz,  and 
fs  =  S\92Hz,  where  fc  is  the  carrier  frequency,  rb  is  the  bit  rate,  and  fs  is  the 

sampling  frequency. 
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Figure  36  Surface  plot  of  the  SCD  estimate  magnitude  for  an  Amplitude  Shift 
Keying  (ASK)  or  a  Binary  Phase  Shift  Keying  (BPSK)  Signal,  using  AUTOFAM, 
with  the  following  parameters:  Af  =  512Hz ,  &a  =  \6Hz ,  fc  =  2048Hz , 
rb=\024Hz,  and  /5  =  8192//z,  where  fc  is  the  carrier  frequency,  rb  is  the  bit 
rate,  and  fs  is  the  sampling  frequency. 
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Figure  37  Contour  plot  of  the  SCD  estimate  magnitude  for  an  Amplitude  Shift 
Keying  (ASK)  or  a  Binary  Phase  Shift  Keying  (BPSK)  Signal,  using  AUTOFAM, 
with  the  following  parameters:  Af  =  5\2Hz,  Aa  =  \6Hz,  fc=204SHz, 
rb=\024Hz,  and  fs  =  8\92Hz,  where  fc  is  the  carrier  frequency,  rb  is  the  bit 
rate,  and  fs  is  the  sampling  frequency. 
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Figure  38  Plots  of  the  SCD  estimate  magnitude  for  an  Amplitude  Shift  Keying 
(ASK)  or  a  Binary  Phase  Shift  Keying  (BPSK)  Signal,  using  AUTOFAM,  for 
a  =  0 ,  a  =  1024Hz ,  a  =  2048/fe ,  a  =  3072Hz  ,  and  a  =  4096Hz ,  respectively,  with 
the  following  parameters:  A/  =  5\2Hz ,  Aa  =  \6Hz ,  fe  =  204SHz ,  rb  =  \024Hz ,  and 
/,  =  8192ife,  where  fc  is  the  carrier  frequency,  rb  is  the  bit  rate,  and  fs  is  the 
sampling  frequency. 
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Figure  39  Surface  plot  of  the  SCD  estimate  magnitude  for  an  Amplitude  Shift 
Keying  (ASK)  or  a  Binary  Phase  Shift  Keying  (BPSK)  Signal,  using  AUTOSSCA, 
with  the  following  parameters:  Af  =  512Hz,  Aa  =  \6Hz,  fc=204SHz, 
rb  =  \024Hz,  and  fs  =  8\92Hz,  where  fe  is  the  carrier  frequency,  rb  is  the  bit 
rate,  and  fs  is  the  sampling  frequency. 
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Figure  40  Contour  plot  of  the  SCD  estimate  magnitude  for  an  Amplitude  Shift 
Keying  (ASK)  or  a  Binary  Phase  Shift  Keying  (BPSK)  Signal,  using  AUTOSSCA, 
with  the  following  parameters:  A/  =  512Hz ,  t\a  =  4  Hz ,  fe  =  204SHz ,  rb  =  1024/fe , 
and  fs  =  8\92Hz ,  where  fc  is  the  carrier  frequency,  rb  is  the  bit  rate,  and  fs  is 
the  sampling  frequency. 
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Figure  41  Plots  of  the  SCD  estimate  magnitude  for  an  Amplitude  Shift  Keying 
(ASK)  or  a  Binary  Phase  Shift  Keying  (BPSK)  Signal,  using  AUTOSSCA,  for 
a  =  0,  a  =  \024Hz,  a  =  204SHz,  a  =  3072Hz,  and  a  =  4096Hz,  respectively,  with 
the  following  parameters:  A/  =  5\2Hz ,  Aa  =  \6Hz ,  fc  =  2048#z ,  rb  =  1024Hz ,  and 
fs  =  S\92Hz,  where  fc  is  the  carrier  frequency,  rb  is  the  bit  rate,  and  /,  is  the 

sampling  frequency. 
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VI.  CONCLUSIONS 

A.  SUMMARY 

The  purpose  of  this  thesis  was  to  implement  the  FAM  and  the  SSCA 
methods  in  MATLAB,  allowing  students,  researchers,  and  engineers  to  take 
advantage  of  the  power  of  the  cyclic  analysis  methods  for  solving  signal 
detection  and  modulation  identification  problems.  Since  MATLAB  is  user  friendly 
and  easily  portable  to  other  operating  systems,  the  implementation  becomes  a 
proving  ground,  easily  modified  and  set  up  on  other  computer  systems. 

The  two  methods  are  implemented  as  user  friendly  as  possible.  The  only 
inputs  required  for  both  are  the  signal(s),  the  sampling  frequency  (fs),  that 
should  be  the  same  for  both,  and  the  desired  resolutions  for  spectrum  frequency 
(/)  and  cyclic  frequency  {a).  The  results  generated  by  both  programs  can  be 
displayed  in  different  ways,  such  as  surface  plots,  contour  plots,  and  cross- 
section  plots. 

Both  programs  generate  a  large  amount  of  data  when  a  good  resolution  is 
desired.  As  a  consequence,  limitations  in  computational  as  well  as  printing 
resources  did  not  allow  the  presentation  of  more  detailed  plots.  One  can  of 
course  easily  create  averaged  spectral  correlation  surfaces,  with  a  coherent  or 
power  averaging  method. 
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B.  SUGGESTIONS 

A  step  that  was  not  possible  to  reach  during  this  work  was  an  analysis 
and  experimentation  of  the  performance  of  cyclic  spectral  analysis  in  a  white- 
noise  contaminated  environment.  A  natural  extension  to  this  work  is  an  analysis 
of  how  well  cyclic  analysis  performs  in  a  white  and  a  colored  noise  environment 
such  as  the  one  imposed  by  the  ocean  to  the  receiving  elements. 

Computationally,  cyclic  spectral  analysis  is  an  expensive  task.  Therefore, 
any  improvements  to  the  speed  of  existing  methods  or  even  completely  new  fast 
algorithms  are  desirable. 

During  the  development  of  this  thesis,  it  was  intended  to  replace  some  of 
the  Fast  Fourier  Transform  (FFT)  operations  by  the  fast  Wavelet  Transform.  A 
follow  on  to  this  work  is  the  replacement  of  the  FFT  operations  by  Wavelet 
Transforms,  and  an  evaluation  of  how  well  the  modified  algorithm  performs  in 
terms  of  identification  and  computational  cost. 
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APPENDIX  A.  CALCULATION  OF  THE  SCD  FUNCTION  OF  AN 
AMPLITUDE-MODULATED  SIGNAL 


Let  us  consider  the  amplitude-modulated  (AM)  signal 

x[n]  =  a[n]  cos\2tu  f  n) ,  (52) 

where  a[n]   is  a  stationary  random  lowpass  signal  with  PSD  Sa(f),  with  no 

spectral  lines  in  its  PSD. 

The  fundamental  parameter  for  second-order  periodicity  in  discrete-time, 
called  cyclic  autocorrelation  function  and  denoted  by   R?[£]  ,  is  given  by 

K¥\=  (x[n]x*[n-  l]e-2Kan)  e'nai .  (53) 

We  can  look  at  R"[e]  from  the  following  point  of  view: 
Let  u[n]  =  x[n]  e-na"  and  v[n]  =  x[n]  e+nran ,  then 

R°[e]=(x[n]e-,7:a"  x'ln-^e'™  e'na() 
=  (x[n\e-*a"  x'[n-£]e-ixa(n-e)) 
=  (x[n]e-'™{x[n-l]e"ra{"-()y 
=  (u[n]v,[n-£]) 

=RAn  (54) 

But,  Since  u[n]  =  x[n]  e~i7ra"  , 
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Ru[£]  =  (u[n]u[n-t]) 

=  [x[n\e-nan  x[n-t]eimln-e)) 
=  (x[«]x*[«-^])  e 


-mat 


e  -incd 


-12*1  ^U 


=  Rx[t]e      K2)  (55) 

and  taking  the  Fourier  transform  leads  to 

Su(f)  =  Sx(f)®s[f  +  ^)=Sx{f  +  ^].  (56) 

Similarly,  if  v[n]  =  x[n]  e'nan ,  then 


i2*|  - 


and 

Sv(f)  =  Sx[f-^j.  (58) 

Thus,  we  can  redefine  the  second-order  periodicity,  by  saying  that  x[n] 
exhibits  second-order  periodicity  if  and  only  if  frequency  translates  (frequency- 
shifted  versions)  of  x[n]  are  correlated  with  each  other. 

As  long  as  the  mean  values  of  the  frequency  translates  u[n]  and  v[n]  are 
zero     (i.  e.,   (w[«])  =  0and   (v[«])  =  0),    x[n]   does  not  contain  finite-amplitude 
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additive  sine-wave  components  at  ±a/2  and  therefore  Sx(f)  has  no  spectral 
lines  at  f  =  ±a/2.  The  crosscorrelation  RmU]  is  actually  a  temporal 
crosscovariance  Kuv[i] ,  that  is, 

KM  =  ({u[n}-(u[n})}{v{n-t]-(v[n-e})Y 

=  (u[n]v*[n-£]) 

=  RM-  (59) 

Based  on  this,  we  can  define  a  temporal  correlation  coefficient,  the 
magnitude  of  which  is  upper  bounded  by  unity,  for  frequency  translates  since 

^  14  =  TTnT  =  7     r  i     r  i\^"  •  (60) 

RM    [kMkM) 

Since  a[«]  is  a  real  stationary  random  (  R"[£]  =  0,  for  all  a  *  0)  signal  with 
zero  mean  (  i.  e.,  (a[n\)  =  0  ),  the  cyclic  autocorrelation  function  of  x[n]  is  given 
by 

R°  [£]  =  (x[n]  x*[n-£]  e''2™ )  e,KCd 

=  (a[n] cos(2/r/>  +  6)  a  [n-  i] co{2tt/0 (n-  i)  +  0]  e"2™ )  eimt 

=  (a[n]a  [n  -  £]- {co^Itt/J)  +  cos[2;r/0 (in  - 1)+  2d]}  e-2™ )  eixai 
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^cos(2tt fj)  {a[n]a[n -  i\e"2^)  e 


mat 


+  -(a[ny[n-  t]cos{27rf0(2n-  £)+20]e-,2*a")e"rae 

i  1/  i[2xfo(2n-e)+20]  -i[2xfa(2»-l)+20] 

= -cos(2x  fj)R:[£]+-  la[n]a-[n-t] 


„-i2xan   \  „mai 

e  e 


=  - Raa  [t\co^27cfj)  +  -ene(a[n\a[n - ^"'■^/>|e'^/.V 


+  —  e      \a\n\a  \n-t\e      y      "'  )  e  y      "' 


=  ±R:[i]co42xfj)  +  ±e™  Rrlf°[tV\e-™  R^[i]. 


(61) 


But  since  a[n]  is  a  stationary  signal  (i.  e.,  R"[£]  =  0  for  all  a  *  0),  the  only 
non-zero  contributions  are  at  a  =  0  and  a  =  ±2f0 .  Thus,  the  cyclic 
autocorrelation  function  becomes 


*?[/]- 


-e»2°Ra[tl  a  =  2f0 


4 


■Ra[£]™(2xfot\       «  =  0 


otherwise 


(62) 


Defining  the  temporal  correlation  coefficient  as 
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sxa[£]= 


(63) 


where  Rx[0]  =  7?;[^]  for  a  =  0,  computed  at  £  =  0 ,  so  #x[0]  =  i?a[0]/2 
So,  S"[i]  becomes 


sxa[£]= 


—  e 
4 


±ne 


R.¥\ 


*M 


=  —e 
2 


±ne 


sM 


a  =  ±2/0 


-Ra[e]M2*fo*) 

-  =  Sa[£]cos(2xfj),       a  =  0. 


\rm 


(64) 


Thus,  the  strength  of  correlation  between  x[n]e~i7rcm  and  x[«-^]  e'"1^"'0  is 


given  by 


l*.M|. 


a  =  ±2f0 


<\sm[£]ux(2xf0t%       a  =  0 

0,  otherwise. 


(65) 


It  can  be  substantial  for  an  amplitude-modulated  signal,  e.g.,|j/[0]|  =  1/2  . 
To  localize  in  the  frequency  domain  the  average  power  (|x[«]|  /  =  i?JO]  in 

a  stationary  signal  x[n],  we  simply  pass  the  signal  x[n]  through  a  set  of 
narrowband  bandpass  filters  and  then  measure  the  average  power  at  the  output 
of  the  filters.  In  the  limit  when  the  bandwidths  of  these  filters  approach  zero,  the 
corresponding   set  of  measurements  of  average   power,    normalized   by  the 
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bandwidth,  constitute  the  "power  spectral  density  (  PSD  )  function",  given  by 

^(y)  =  lim -^(^[/i] ®  x[«]  |2) ,  (66) 

where  hfB[n]  is  the  discrete-impulse  response  of  a  bandpass  filter  with  center 
frequency  / ,  bandwidth  B ,  and  unity  gain  at  the  band  center. 

In  a  similar  fashion,  to  localize  in  the  frequency  domain  the  correlation 

(w[«]v*[«])  =  (|x[«]|  e~'2'ran)^  R"[0]  of  frequency-shifted  versions  u[n]  and  v[n]  of 

a  cyclostationary  signal  x[n] ,  we  pass  both  of  the  two  frequency  translates  u[n] 
and  v[n]  of  x[n]  through  the  same  set  of  bandpass  filters.  Again,  in  the  limit 
when  the  bandwidth  of  these  filters  approach  zero,  the  corresponding  set  of 
measurements  of  the  temporal  correlation  of  the  filtered  signals  constitute  the  so 
called  "spectral-correlation  density  (SCD  )  function",  given  by 

5"(/")  sjJ2^(feW®  «[»]}{**  W®  *]}*)•  <67> 

that  yields  the  spectral  density  of  correlation  between  u[n]  and  v[n]  at  frequency 
/ ,  which  is  identical  to  the  spectral  density  of  correlation  in  x[n]  at  frequencies 
/  +  a/2  and  f-a/2. 

It's  a  well  known  result,  called  "Wiener  relation"  as  opposed  to  its 
probabilistic  counterpart  called  "Wiener-Khinchin  relation",  that  the  PSD  is  equal 
to  the  Fourier  transform  of  the  autocorrelation  function, 
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«,(/)=  I  A.M^-  (68) 


f=-oo 


Similarly,  the  SCD  can  be  obtained  by  Fourier  transforming  the  "cyclic 
autocorrelation  function", 

*;(/)=  E*;M*"'2jr/r-  (69) 

^=-00 

This  result  is  known  as  the  "cyclic  (  periodic  )  Wiener  relation". 

We  observe  that  since   R"U]   is  periodic  in   a  with  period  two,  so  is 

K+2  M  =  (*M *'[n-  £]e-i2*{a+2)" )  e'*{a+2)e 

=  (x[n]xt[n-e]e-i2*a"  e-4*")  e'nai  e'l7d  ,  (70) 

and,  since  n  and  £  are  integers, 

«-/4-=«/2*€=ll  (71) 

so, 

Rr2[^={x[n]x'[n-e]e-l2^)e"a( 

=  RaM\.  (72) 

Thus, 
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sr2(/)=X  Rr2[*]e-i2*j 


Tft 

(=-0 


=  £  R?[t]e-i2*fi 

=  S;(f).  (73) 

Also,  since  I  takes  on  only  integer  values,  then  S°(f)  is  periodic  in  / 
with  period  one 

S,°(/  +  l)=E  R:[£]e-i2*f+l)t 

=  Z  KV\*~n*Si  e'1™  -  (74) 

and,  since  ^  is  an  integer,  we  have 

e~i2!re=\,  (75) 

then, 

=  5:(/).  (76) 

Furthermore,  S°(f)  also  exhibits  the  following  periodicity: 

Srl{f  +  V2)  =  S;(f).  (77) 
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This  is  easily  shown  by 


sr'|/+i)  =  £*r'M*-'2*(/4> 


=  |;{*["|x'[»-<k'2*(a*,,")<;'*l°+')'  e-1""  e 


■iid 


(=-k 


=  Yi(x[n]x'[n-l]e-i2m  e-i2m)eixae  eM  e-2nfi  e 


-inl 


=  Yd{x[n]x'[n-i\e-am  e~i2m)  einat  e~aKfl  , 


£=-CO 


and,  since  n  is  an  integer,  we  have 


(78) 


e"2™  =  1 , 


(79) 


hence, 


=  I  KW-*** 


(=-<*> 


=  saAf) 


(80) 


We  also  define  a  spectral  correlation  coefficient,  p"(f) ,  given  by 


p!{/)- 


s;(f) 


*j/+g</-t 
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(81) 


Since  |>o"(/)|  is  bounded  by  [o,l],  it  is  a  convenient    measure  of    the 

degree  of  local  spectral  redundancy  that  results  from  spectral  correlation. 

Going    back   to   the   AM    signal,    by    Fourier   transforming    the    cyclic 
autocorrelation  function,  we  get 


s;{/)= 


\e^%(/),  a  =  ±2f„ 

7s.{f+f.)+7s-(f-f'\       a  =  0 


otherwise 


(82) 


and  so,  the  spectral  correlation  coefficient  is  given  by: 


H/+f J =  ^W-for  a  =  °'  comPuted  at  /+f 


=  jSa{f  +  fo)+jSa{f-f0),  computed  at  f+% 


1      (        a         )      1  a 

^Sa[f  +  -  +  f0J+-Sa[f  +  --f0 


(83) 


Also, 


a)      1      (        a        \     1      (        a 


(84) 


So  that,  fora  =  2/0, 
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s{f  +  ^=Sx{f  +  f0)  =  \sa{f  +  2f0)  +  \sa{f), 


(85) 


and, 


s{f  ~^=sAf  -  L)  =  \sXf)  +  \sXf  -2f0), 


and  for  a  =  -2f0 , 


and, 


and,  for  a  =  0, 


and, 


S,(/  +  f)=S,(/-/„)  =  is„(/)  +  ^(/-2/„), 


/-f)=S,(/  +  /„)  =  isa(/  +  2/,)+|s„(/), 


S,[/  +  f)  =  «,(/)  =  jSo(/  +  /„)+is„(/-/0), 


^(/-f)  =  ^(/)  =  ^s„(/+/„)+^s„(/-/„), 


then, 


fW- 


s.(fy™ 


f,  a  =  ±2/, 

{[s„  (/  +  2/0 )  +  5a  (/)  Js,  (/)  +  Sa{f  -  2/0 )]} 2 

1,  a  =  0 

0,  otherwise. 


(86) 


(87) 


(88) 


(89) 


(90) 


(91) 
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Thus,  the  strength  of  correlation  between  the  spectral  components  of  x[n] 
at  frequencies  f  +  a/2  and  f-a/2  is  unity  (i.  e.,  |p"(/)|  =  l,  for  \f\<f0  and 
a  =  ±2f0),  provided  that  a[n]   is  bandlimited  to   \f\<f0.  Hence   Sa(f)  =  0  for 

\Mfo- 
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APPENDIX  B.  FUNCTION  AUTOFAM 


function  [Sx, alphao, fo] =autofam(x, f s, df , dalpha) 

%  AUTOFAM (X,FS,DF, DALPHA)  computes  the  spectral  auto- correlation 

%  density  function  estimate  of  the  signal  X,  by  using  the  FFT 

%  Accumulation  Method (FAM) .  Make  sure  that  DF  is  much  bigger 

%  than  DALPHA  in  order  to  have  a  reliable  estimate. 

% 

%  INPUTS : 

%  X       -  input  column  vector; 

%  FS     -  sampling  rate; 

%  DF     -  desired  frequency  resolution;  and 

%  DALPHA  -  desired  cyclic  frequency  resolution. 

Q, 
O 

%  OUTPUTS : 

%  SX     -  spectral  correlation  density  function  estimate; 

%  ALPHAO  -  cyclic  frequency;  and 

%  FO     -  spectrum  frequency. 

Q. 
O 

%  Author:  E.L.Da  Costa, 9/28/95 . 

if  nargin  ~=  4 

error (v Wrong  number  of  arguments.'); 
end 


%  Definition  of  Parameters 


Np=pow2 (nextpow2 (fs/df ) ) ;      %  Number  of  input  channels,  defined 

%  by  the  desired  frequency 
%  resolution (df)  as  follows: 
%  Np=fs/df,  where  fs  is  the  original 
%  data  sampling  rate.  It  must  be  a 
%  power  of  2  to  avoid  truncation  or 
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L=Np/4 ; 


P=pow2 (nextpow2 (f s/dalpha/L) 


N=P*L; 


%  zero-padding  in  the  FFT  routines; 

%  Offset  between  points  in  the  same 

%  column  at  consecutive  rows  in  the 

%  same  channelization  matrix.  It 

%  should  be  chosen  to  be  less  than 

%  or  equal  to  Np/4 ; 

%  Number  of  rows  formed  in  the 

%  channelization  matrix,  defined 

%  by  the  desired  cyclic  frequency 

%  resolution (dalpha)  as  follows: 

%  P=f s/dalpha/L. It  must  be  a  power 

%  of  2; 

%  Total  number  of  points  in  the 

%  input  data . 


%  Input  Channelization  % 

oooooooooooooooooooooooo 

if  length (x) <N 

x(N)=0; 
elseif  length (x) >N 

x=x(l:N) ; 
end 

NN= (P-l) *L+Np; 
XX=X; 
Xx(NN)=0; 
xx=xx ( : ) ; 
X=zeros (Np, P) ; 
for  k=0:P-l 

X( : ,k+l) =xx(k*L+l:k*L+Np) ; 
end 


Windowing 
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a=hamming (Np) ; 
XW=diag(a) *X; 
XW=X  ; 


%  First  FFT  % 

2-2-9-9-2-S-S-9-2-2-9-9-9- 
ooooooooooooo 

XFl=f f t (XW) ; 

XFl=fftshift (XF1) ; 

XF1= [XF1( : ,P/2+l:P)  XF1 ( : ,l:P/2) ] ; 


Downconver s  ion 


E= zeros (Np, P) ; 
for  k=-Np/2:Np/2-l 

for  m=0:P-l 

E(k+Np/2+l,m+l) =exp(-i*2*pi*k*m*L/Np) ; 

end 
end 

XD=XF1.*E; 
XD=conj (XD' ) ; 


%  Multiplication  % 

000000000000000*00*0 

XM=zeros  (P,Np'K2)  ; 
for  k=l:Np 

for  1=1 :Np 

XM(: , (k-l)*Np+l)=(XD(: ,k) .*COnj (XD ( : ,1) ) ) ; 

end 
end 
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%  Second  FFT  % 

2-S-S-B-B-2-S-2-S-9-S-S-2-2- 

oooooooooooooo 

XF2=fft (XM) ; 

XF2=fftshift (XF2) ; 

XF2= [XF2 ( : , Np^2/2+l :Np^2 )     XF2 ( : ,l:NpA2/2) ] ; 

XF2=XF2 (P/4 : 3*P/4 , : ) ; 

M=abs (XF2) ; 

alphao=-l:l/N:l; 

fo=-.5:l/Np: .5; 

Sx=zeros (Np+1, 2*N+1) ; 

for  kl=l: P/2+1 

for  k2=l:Np*2 

if  rem(k2,Np) ==0 

l=Np/2-l; 
else 

l=rem(k2,Np) -Np/2-1; 
end 

k=ceil (k2/Np) -Np/2-1; 
p=kl-P/4-l; 

alpha= (k-1) /Np+ (p-1) /L/P; 
f=(k+l)/2/Np; 
if  alpha<-l    alpha>l 

k2=k2+l; 
elseif  f <- . 5    f > . 5 

k2=k2+l; 
else 

kk=l+Np* (f+.5) ; 
11=1+N* (alpha+1) ; 
Sx(kk,ll) =M(kl,k2) ; 
end 
end 
end 
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APPENDIX  C.   FUNCTION  AUTOSSCA 


function  [Sx, alphao, fo] =autossca (x, f s, df , dalpha) 

%  AUTOSSCA (X, FS,DF, DALPHA)  computes  the  spectral  auto- 

%  correlation  density  function  estimate  of  the  signal  X, 

%  by  using  the  Strip  Spectral  Correlation  Algorithm  (SSCA) 

%  Make  sure  that  DF  is  much  bigger  than  DALPHA  in  order  to 

%  have  a  reliable  estimate. 

% 


INPUTS : 

X 

FS 

DF 

DALPHA 


input  column  vector; 

sampling  rate; 

desired  frequency  resolution;  and 

desired  cyclic  frequency  resolution. 


OUTPUTS : 

SX     -  spectral  auto-correlation  density  function  estimate; 

ALPHAO  -  cyclic  frequency;  and 

FO     -  spectrum  frequency. 

Author:  E.L.Da  Costa, 9/28/95 . 


if  nargin  ~=  4 

error (* Wrong  number  of  arguments'); 
end 


Definition  of  Parameters 


Np=pow2 (nextpow2 (fs/df ) ) ; 


%  Number  of  input  channels,  defined 

%  by  the  desired  frequency 

%  resolution (df)  as  follows: 

%  Np=fs/df,  where  fs  is  the  original 

%  data  sampling  rate.  It  must  be  a 
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L=Np/4 ; 


P=pow2 (nextpow2 (fs/dalpha/L) ) ; 


N=P*L; 


power  of  2  to  avoid  truncation  or 
zero-padding  in  the  FFT  routines; 
Offset  between  points  in  the  same 
column  at  consecutive  rows  in  the 
same  channelization  matrix.  It 
should  be  chosen  to  be  less  than 
or  equal  to  Np/4; 
Number  of  rows  formed  in  the 
channelization  matrix,  defined  by 
the  desired  cyclic  frequency 
resolution (dalpha)  as  follows: 
P=f s/dalpha/L.  It  must  be  a  power 
of  2; 

Total  number  of  points  in  the 
input  data . 


%  Input  Channelization  % 

9-9-9-9-9^S-2-9-9-9-9-2-2-S-2-2-2-2-2-9-S-B-S-2- 
"00000*0000000000000000000 

if  length (x) <N 

x(N)=0; 

disp('you  will  not  get  the  desired  resolution  in  cyclic  frequency'); 

dalpha=f s/N; 

disp (  [ ' cyclic  frequency  resolution= ■ , num2str (dalpha) ] )  ; 
elseif  length (x) >N 

x=x(l:N) ; 
end 

NN= (P-l) *L+Np; 
xx=x  ; 
XX (NN) =0; 
xx=xx ( : ) ; 
X=zeros (Np, P) ; 
for   k=0:P-l 

X( : ,k+l) =xx(k*L+l:k*L+Np) ; 
end 
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%  Windowing  % 

S-9-2-9-9-9-9-9-2-2-S-S-S- 

OOOOOOOOOOOOO 

a=hamming (Np) ; 
XW=diag(a) *X; 

S-S-S-S-S-S-S-S-S-S-S-S-S- 
000000000*0000 

%  First  FFT  % 

ooooooooooooo 

XFl=fft (XW) ; 

XFl=fftshift (XF1) ; 

XF1=[XF1 ( : ,P/2+l:P)  XF1 ( : ,l:P/2)] ; 


%  Downconversion 


E= zeros (Np, P) ; 

for  k=-Np/2:Np/2-l 

for  m=0:P-l 

E(k+Np/2+l,m+l)=exp(-i*2*pi*k*m*L/Np) ; 

end 
end 
XD=XF1.*E; 


%  Replication  % 

9-9-9-S-9-S-9-S-2-9-2-?—  2-S-S- 
ooooooooooooooo 

XR=zeros (Np, P*L) ; 
for  k=l:P 

XR( : , (k-1) *L+l:k*L)=XD(:,k)*ones(l,L) ; 

end 


%  Multiplication  % 
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xc=ones (Np, 1) *x' ; 
XM=XR . *xc ; 
XM=conj (XM1  ) ; 


Second  FFT 


XF2=fft (XM) ; 

XF2=fftshift (XF2) ; 

XF2=[XF2 ( : ,Np/2+l:Np)  XF2 ( : ,l:Np/2) ] ; 

M=abs(XF2) ; 

alphao= (-1:1/N:1) *fs; 

fo=(- .5:1/Np: .5) *fs; 

Sx=zeros (Np+1, 2*N+1) ; 

for  kl=l:N 

for  k2=l:Np 

alpha= (kl-1) /N+ (k2-l) /Np-1; 

f=((k2-l)/Np- (kl-l)/N)/2; 

k=l+Np* (f+. 5) ; 

1=1+N* (alpha+1) ; 

Sx(k,l)=M(kl,k2) ; 
end 
end 
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APPENDIX  D.   FUNCTION  CROSSFAM 


function  [Sxy, alphao, f o] =crossf am (x, y, f s , df , dalpha) 

%        CROSSFAM (X, Y, FS,DF, DALPHA)  computes  the  spectral  cross- 

%        correlation  density  function  estimate  of  the  signals  X 

%        and  Y,  by  using  the  FFT  Accumulation  Method (FAM) .  Make 

%        sure  that  DF  is  much  bigger  than  DALPHA  in  order  to  have 

%        a  reliable  estimate. 

% 

% 

% 

% 

% 

% 

% 

Q, 
O 

% 

Q, 

O 

% 

"o 

a, 
O 


INPUTS : 

X 

Y 

FS 

DF 

DALPHA 


input  column  vector; 

input  column  vector; 

sampling  rate; 

desired  frequency  resolution;  and 

desired  cyclic  frequency  resolution. 


OUTPUTS : 

SXY    -  spectral  cross-correlation  density  function  estimate; 

ALPHAO  -  cyclic  frequency;  and 

FO      -  spectrum  frequency. 

Author:  E.L.Da  Costa, 9/28/95 . 


If  nargin  ~=  5 

error ( 'Wrong  number  of  arguments . ' ) ; 
end 


%  Definition  of  Parameters 


Np=pow2 (nextpow2 (fs/df ) ) ; 


%  Number  of  input  channels,  defined 

%  by  the  desired  frequency 

%  resolution (df)  as  follows: 

%  Np=fs/df,  where  fs  is  the  original 
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L=Np/4 ; 


P=pow2 (nextpow2 (fs/dalpha/L) ) ; 


N=P*L; 


data  sampling  rate.  It  must  be  a 
power  of  2  to  avoid  truncation  or 
zero-padding  in  the  FFT  routines ; 
Offset  between  points  in  the  same 
column  at  consecutive  rows  in  the 
same  channelization  matrix.  It 
should  be  chosen  to  be  less  than 
or  equal  to  Np/4; 
Number  of  rows  formed  in  the 
channelization  matrix,  defined  by 
the  desired  cyclic  frequency 
resolution (dalpha)  as  follows: 
P=f s/dalpha/L.  It  must  be  a  power 
of  2; 

Total  number  of  points  in  the 
input  data . 


%  Input  Channelization  % 

9-9-9-S-9-9-S-B-9-S-S-S-9-2-2-2-B-S-B-9-2-S-2-9- 
oooooooooooooooooooooooo 

if  length (x) <N 

x(N)=0; 
elseif  length (x) >N 

x=x(l:N) ; 
end 
if  length (y) <N 

y(N)=0; 
elseif  length (y) >N 

y=y(l:N) ; 
end 

NN= (P-l) *L+Np; 
XX=x; 

yy=y; 

xx (NN) =0; 
yy(NN)=0; 


104 


xx=xx ( : ) ; 

yy=yy ( : ) ; 
X=zeros (Np, P) ; 
Y= zeros (Np, P) ; 
for  k=0:P-l 

X( : ,k+l) =xx(k*L+l:k*L+Np) ; 

Y( : , k+l)=yy(k*L+l:k*L+Np) ; 
end 


%  Windowing  % 

S-S-S-S-S-S-S-S-S-S-S-S-9- 

ooooooooooooo 

a=hamming (Np) , 
XW=diag(a) *X; 
YW=diag(a) *Y; 


%  First  FFT  % 

ooooooooooooo 

XFl=f f t (XW) ; 

YFl=f ft  (YW)  ; 

XFl=fftshift (XF1) ; 

YFl=fftshift (YF1) ; 

XF1=[XF1( : ,P/2  +  l:P)  XF1 ( : ,l:P/2)  ]  ; 

YF1=[YF1 ( : ,P/2  +  l:P)  YF1 ( : ,l:P/2)  ]  ; 


%  Downconversion 


E=zeros (Np, P) ; 

for  k=-Np/2+l:Np/2 

for  m=0:P-l 

E(k+Np/2,m+l) =exp(-i*2*pi*k*m*L/Np) ; 

end 
end 
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XD=XF1.*E; 
YD=YF1.*E; 
XD=COnj (XD ' ) ; 
YD=conj (YD1 ) ; 


%  Multiplication  % 

2-9-S-9-2-9-9-S-2-9-2-S-S-2-2-S-B-S- 
oooooooooooooooooo 

XYM=zeros (P,NpA2) ; 
for  k=l:Np 

for  1=1 :Np 

XYM( : , (k-1) *Np+l) =XD ( : ,k) . *conj (YD ( : ,1) ) ; 

end 
end 


%  Second  FFT 


XYF2=fft (XYM) ; 

XYF2=fftshift (XYF2) ; 

XYF2=[XYF2 (: , Np^2/2+l :Np"2)  XYF2 (: ,l:Np*2/2) ] ; 

XYF2=XYF2 (P/4:3*P/4, :) ; 

M=abs(XYF2) ; 

alphao= ( - 1 : 1 /N : 1 ) *  f s ; 

fo= (- .5:1/Np: .5) *fs; 

Sxy=zeros (Np+1, 2*N+1) ; 

for  kl=l:P/2+l 

for  k2=l:NpA2 

if  rem(k2,Np) ==0 
l=Np/2; 

else 

l=rem(k2,Np) -Np/2 ; 

end 

k=ceil (k2/Np) -Np/2; 

106 


p=kl-P/4-l; 

alpha= (k-1) /Np+ (p-1) /L/P; 

f=(k+l)/2/Np; 

if  alpha<-l    alpha>l 

k2=k2+l; 
elseif  f<-.5  |  f>.5 

k2=k2+l; 
else 

kk=l+Np* (f+.5) ; 

11=1+N* (alpha+1) ; 

Sxy(kk,ll)=M(kl,k2) ; 
end 


end 


end 
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APPENDIX  E.   FUNCTION  CROSSSSCA 


function 
% 

9- 
o 


% 

to  havi 

% 

% 

INPUTS 

% 

X 

% 

Y 

% 

FS 

o 

DF 

% 

DALPHA 

[Sxy, alphao, f o] =crossssca (x, y, f s , df , dalpha) 
CROSSSSCA (X,Y,FS,DF, DALPHA)  computes  the  spectral  cross- 
correlation  density  function  estimate  of  the  signals  X 
and  Y,  by  using  the  Strip  Spectral  Correlation  Algorithm 
(SSCA) .  Make  sure  that  DF  is  much  bigger  than  DALPHA  in  order 
to  have  a  reliable  estimate. 


input  column  vector; 

input  column  vector; 

sampling  rate; 

desired  frequency  resolution;  and 

desired  cyclic  frequency  resolution. 


%  OUTPUTS : 

%  SXY     -  spectral  cross-correlation  density  function  estimate; 

%  ALPHAO  -  cyclic  frequency;  and 

%  FO     -  spectrum  frequency. 

O 

%        Author:  E.L.Da  Costa, 9/28/95 . 

If  nargin  -=  5 

error ( 'Wrong  number  of  arguments . ' ) ; 
end 


Definition  of  Parameters 


Np=pow2 (nextpow2 (fs/df ) ) ; 


%  Number  of  input  channels,  defined 

%  by  the  desired  frequency 

%  resolution (df)  as  follows: 

%  Np=fs/df,  where  fs  is  the  original 
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L=Np/4 ; 


P=pow2 (nextpow2 (fs/dalpha/L) ) ; 


N=P*L; 


data  sampling  rate.  It  must  be  a 
power  of  2  to  avoid  truncation  or 
zero-padding  in  the  FFT  routines; 
Offset  between  points  in  the  same 
column  at  consecutive  rows  in  the 
same  channelization  matrix.  It 
should  be  chosen  to  be  less  than 
or  equal  to  Np/4; 
Number  of  rows  formed  in  the 
channelization  matrix,  defined  by 
the  desired  cyclic  frequency 
resolution (dalpha)  as  follows: 
P=f s/dalpha/L.  It  must  be  a  power 
of  2; 

Total  number  of  points  in  the 
input  data . 


%  Input  Channelization  % 

9-9-2-S-S-S-9-S-S-9-2-S-S-2-9-2-S-2-S-S-S-S-S-9- 
oooooooooooooooooooooooo 

if  length (x) <N 

x(N)=0; 
elseif  length (x) >N 

x=x(l:N) ; 
end 
if  length (y) <N 

y(N)=0; 
elseif  length (y) >N 

y=y(l:N) ; 
end 

NN=(P-1) *L+Np; 
xx=X; 
Xx(NN)=0; 
xx=xx ( : ) ; 
X=zeros (Np, P) ; 


110 


for  k=0:P-l 

X( : ,k+l)=xx(k*L+l:k*L+Np) ; 
end 


9-9-S-9-2-2-9-S-S-S-9-2-9- 

ooooooooooooo 


%  Windowing 


9-2-S-2-S-9-9-9-2-S-S-9-2- 
ooooooooooooo 


a=hamming (Np) ; 
XW=diag(a) *X; 


%  First  FFT 


XFl=fft (XW) ; 

XFl=fftshift (XF1) ; 

XF1=[XF1 ( : ,P/2+l:P)  XF1 ( :,l:P/2)] ; 


Downconversion  % 


E=zeros (Np, P) ; 
for  k=-Np/2+l:Np/2 

for  m=0:P-l 

E(k+Np/2,m+l) =exp(-i*2*pi*k*m*L/Np) ; 

end 
end 

XD=XF1.*E; 
XD=conj (XD1 ) ; 


%  Replication  % 


ooooooooooooooo 


XR=zeros (Np,P*L) ; 
for  k=l:P 
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XR( : , (k-1) *L+l:k*L) =XD ( : ,k) *ones (1,L) ; 
end 


%  Multiplication  % 


ooooooooooooooo 


yc=ones (Np, 1) *y ' ; 
XYM=XR . *yc ; 
XYM=COnj (XYM' ) ; 


%  Second  FFT 


%  Second  FFT  % 

S-S-S-2-2-S-9-S-9-2-S-2-2-9- 
oooooooooooooo 

XYF2=fft (XYM) ; 

XYF2=fftshift (XYF2) ; 

XYF2=[XYF2 ( : , Np^2/2+l :Np"2 )  XYF2 ( : ,l:NpA2/2) ] ; 

M=abs (XYF2) ; 

alphao= (-1 :1/N:1) *fs; 

fo=(-.5:l/Np: .5)*fs; 

Sxy=zeros (Np+1, 2*N+1) ; 

for  kl=l:N 

for  k2=l:Np 

alpha= (kl-i) /N+ (k2-l) /Np-l; 

f=( (k2-l)/Np-(kl-l)/N)/2; 
k=l+Np* (f+.5) ; 
1=1+N* (alpha+1) ; 
Sxy(k, 1) =M(kl,k2) ; 
end 
end 
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APPENDIX  F.   PLOTTING  ROUTINES 


rSurface    Plot5 


figure (1) 

surf 1 (alphao, fo, Sx) ; 

view(-37.5,60) ; 

title ('SCD  estimate  using  FAM'); 

xlabel ( ' alpha ' ) ; 

ylabel  (  ' f ' )  ; 

zlabel ( 'Sx' ) ; 


%Contour  Plot% 


figure (2) 

contour (alphao, fo, Sx) ; 
xlabel ( *  alpha (Hz) ' ) ; 
ylabel ('f (Hz) ' ) ; 


%Cross-Section  Plots% 


figure (3) 

plot (fo, Sx ( : , l+N* (alpha/f s+1) ) ) ;    %  alpha  is  the  desired  cyclic 

%  frequency, 
xlabel ( %f (Hz) ' ) ; 
ylabel ( x  Sx (alpha) ' ) ; 
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